Project 0.1

#define _CRT_SECURE_NO_DEPRECATE
#include <stdio.h>
#include <conio.h>
#include <iostream>
#include <math.h>
#include <time.h>
#include <fftw3.h>

#define M_PI 3.1415926535897932384626433832795
double alph = 0.4;

double t = 0.001;

using namespace std;

void main (void) {
	int i, M = 7;
    const int dimArr = int(pow(2.0, M));
	double dx = 1.0/dimArr;
	double xCur;

	FILE *ini, *res, *fk_data, *Euler;
	ini = fopen("initial.dat", "w");
	fk_data = fopen("fk.dat", "w");
	res = fopen("result.dat", "w");
	Euler = fopen("Euler.txt", "w");
	
	fftw_complex *fx, *fk;	//Для вычисления линейных членов
	fftw_plan fx_2_fk, fk_2_fx_;
	fx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	fk = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	fx_2_fk = fftw_plan_dft_1d(dimArr, fx, fk, FFTW_FORWARD, FFTW_ESTIMATE);
	fk_2_fx_ = fftw_plan_dft_1d(dimArr, fk, fx, FFTW_BACKWARD, FFTW_ESTIMATE);
	
	fftw_complex *fx2, *fk2;	//Для вычисления нелинейных членов
	fftw_plan fx2_2_fk2, fk2_2_fx2_;
	fx2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	fk2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	fx2_2_fk2 = fftw_plan_dft_1d(dimArr, fx2, fk2, FFTW_FORWARD, FFTW_ESTIMATE);
	fk2_2_fx2_ = fftw_plan_dft_1d(dimArr, fk2, fx2, FFTW_BACKWARD, FFTW_ESTIMATE);
	
	fftw_complex *k1x, *k1k;	//Для вычисления K1 из Рунге-Кутты
	fftw_plan k1x_2_k1k, k1k_2_k1x_;
	k1x = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	k1k = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	k1x_2_k1k = fftw_plan_dft_1d(dimArr, k1x, k1k, FFTW_FORWARD, FFTW_ESTIMATE);
	k1k_2_k1x_ = fftw_plan_dft_1d(dimArr, k1k, k1x, FFTW_BACKWARD, FFTW_ESTIMATE);
	

	for (i=0; i<dimArr; i++)		//Задание массива линейных членов
	{
		xCur = 2.0*M_PI*i*dx;
		fx[i][0] = 0.1*sin(xCur);
		fx[i][1] = 0.0;
		fprintf(ini,"%1.5f   %1.5f   %1.5f\n",xCur, fx[i][0],fx[i][1]);
	}

	fftw_execute(fx_2_fk);	//Прямое преобразование линейных членов
	for (i=0; i<dimArr; i++){
		fk[i][0] = fk[i][0]/(double)dimArr;
		fk[i][1] = fk[i][1]/(double)dimArr;
		fprintf(fk_data,"%d   %1.5f   %1.5f\n",i, fk[i][0],fk[i][1]);
	}

	for (i=0; i<dimArr; i++)		//Задание массива нелинейных членов
	{
		xCur = 2.0*M_PI*i*dx;
		fx2[i][0] = 0.01*sin(xCur)*sin(xCur); //U^2
		fx2[i][1] = 0.0;
	}

	fftw_execute(fx2_2_fk2);	//Прямое преобразование нелинейных членов
	for (i=0; i<dimArr; i++){
		fk2[i][0] = fk2[i][0]/(double)dimArr;
		fk2[i][1] = fk2[i][1]/(double)dimArr;
	}

	//Cчитаем производные линейные в пространстве Фурье
	fk[0][0] = 0.0;		
	fk[0][1] = 0.0;
	fk[dimArr/2][0] = 0.0;
	fk[dimArr/2][1] = 0.0;

	for (i=1; i<dimArr/2; i++)		
	{
		xCur = 2.0*M_PI*i*dx;
		fk[i][0] = fk[i][0] * (alph*alph*i*i - alph*alph*alph*alph*i*i);
		fk[i][1] = fk[i][1] * (alph*alph*i*i - alph*alph*alph*alph*i*i);

		fk[dimArr - i][0] = fk[dimArr - i][0] * (alph*alph*i*i - alph*alph*alph*alph*i*i);
		fk[dimArr - i][1] = fk[dimArr - i][1] * (alph*alph*i*i - alph*alph*alph*alph*i*i);
	}

	//Считаем производные нелинейные в пространстве Фурье
	fk2[0][0] = 0.0;		
	fk2[0][1] = 0.0;
	fk2[dimArr/2][0] = 0.0;
	fk2[dimArr/2][1] = 0.0;

	for (i=1; i<dimArr/2; i++)		
	{
		fk2[i][0] = fx[i][1] * (alph*i/2);
		fk2[i][1] = fx[i][0] * (-alph*i/2);

		fk2[dimArr - i][0] = fx[dimArr - i][0] * (-alph*i/2);
		fk2[dimArr - i][1] = fx[dimArr - i][1] * (alph*i/2);
	}
	
	//Считаем К1 для Рунге-Кутты в пространстве Фурье
	for (i=0; i<dimArr; i++)		
	{
		k1k[i][0] = (fk[i][0] + fk2[i][0])*t;
		k1k[i][1] = (fk[i][1] + fk2[i][1])*t;
	}
	//Переходим обратно в пространство координат для K1
	fftw_execute(k1k_2_k1x_);

	//попробуем просто U1 = U0 + k1 (Эйлер)

	for (i=0; i<dimArr; i++)		
	{
		xCur = 2.0*M_PI*i*dx;
		fprintf(Euler, "%1.5f   %1.5f	%1.5f\n",xCur,k1x[i][0] + 0.1 * sin(xCur), k1x[i][1]);
	}
	
	
	
	
	
	
	
	
	
	
	fftw_destroy_plan(fx_2_fk);
	fftw_destroy_plan(fk_2_fx_);
	fftw_free(fx);
	fftw_free(fk);

	fftw_destroy_plan(fx2_2_fk2);
	fftw_destroy_plan(fk2_2_fx2_);
	fftw_free(fx2);
	fftw_free(fk2);
	
	fftw_destroy_plan(k1x_2_k1k);
	fftw_destroy_plan(k1k_2_k1x_);
	fftw_free(k1x);
	fftw_free(k1k);
	
	fclose(ini);
	fclose(res);
	fclose(fk_data);
	fclose(Euler);

	getchar();
}



/*fftw_execute(fk_2_fx_); //Обратное преобразование
	for (i=0; i<dimArr; i++){
		xCur = 2.0*M_PI*i*dx;
		fprintf(res,"%1.5f   %1.5f   %1.5f\n", xCur, fx[i][0], fx[i][1]);
	}*/