Project_FINAL

gistfile1.cpp

#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.45;

double t = 0.001;
double TIME = 180.;
int M = 4;
const int dimArr = int(pow(2.0, M));

using namespace std;

void f_nonlinear(fftw_complex * u, fftw_complex * k, int dimArr,fftw_complex * fk2)
{
	double value;
	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);

	for (int i=0; i<dimArr; i++)		//Сложили U + Ki*t/2
	{
		if (k == NULL)
		{

			fk2[i][0] = u[i][0];
			fk2[i][1] = u[i][1];
		}	
		else{

			fk2[i][0] = u[i][0] + t*k[i][0]/2;
			fk2[i][1] = u[i][1] + t*k[i][1]/2;}
	}


	fftw_execute(fk2_2_fx2_); 
	for (int i=0; i<dimArr; i++)		// ^2
	{
		value = fx2[i][0];
		fx2[i][0] = fx2[i][0]*fx2[i][0] - fx2[i][1]*fx2[i][1];
		fx2[i][1] = 2*value*fx2[i][1];
	}

	fftw_execute(fx2_2_fk2);	//Прямое преобразование f2x -> f2k
	for (int i=0; i<dimArr; i++){
		fk2[i][0] = fk2[i][0]/(double)dimArr/2;
		fk2[i][1] = fk2[i][1]/(double)dimArr/2;
	}

	//Пересчет производных
	fk2[0][0] = 0.0;		
	fk2[0][1] = 0.0;
	fk2[dimArr/2][0] = 0.0;
	fk2[dimArr/2][1] = 0.0;

	for (int i=1; i<dimArr/2; i++)
	{
		if (k == NULL)
		{
			value = fk2[i][0];
			fk2[i][0] = (u[i][0])*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) + fk2[i][1] * i * alph;// 
			fk2[i][1] = (u[i][1])*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) - value * i * alph;//

			value = fk2[dimArr - i][0];
			fk2[dimArr - i][0] = (u[dimArr - i][0])*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) + fk2[dimArr - i][1] * (-i) * alph;// 
			fk2[dimArr - i][1] = (u[dimArr - i][1])*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) - value * (-i) * alph;// 

		}		
		else{
			value = fk2[i][0];
			fk2[i][0] =(u[i][0] + t*k[i][0]/2)*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) + fk2[i][1] * i * alph;//  
			fk2[i][1] = (u[i][1] + t*k[i][1]/2)*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) - value * i * alph; //

			value = fk2[dimArr - i][0];
			fk2[dimArr - i][0] = (u[dimArr - i][0] + t*k[dimArr - i][0]/2)*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) +  fk2[dimArr - i][1] * (-i) * alph; //
			fk2[dimArr - i][1] = (u[dimArr - i][1] + t*k[dimArr - i][1]/2)*(pow(alph,2)*pow((double)i,2)-pow(alph,4)*pow((double)i,4)) - value * (-i) * alph;//
		}
	}

	fftw_destroy_plan(fx2_2_fk2);
	fftw_destroy_plan(fk2_2_fx2_);
	fftw_free(fx2);
}
void Aliasing(fftw_complex *temp){
	/*for(int i=dimArr/4; i<dimArr*3/4; i++){
		temp[i][0] = 0.0;
		temp[i][1] = 0.0;
	}*/
}

void main (void) {
	int i;
	double dx = 1.0/dimArr;
	double xCur;

	double tCur = 0.0;

	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 *k1x, *k1k;	//Для вычисления K1 из Рунге-Кутты
	k1x = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	k1k = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);

	fftw_complex *k2x, *k2k;	//Для вычисления k2 из Рунге-Кутты
	k2x = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	k2k = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);

	fftw_complex *k3x, *k3k;	//Для вычисления k3 из Рунге-Кутты
	k3x = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	k3k = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);

	fftw_complex *k4x, *k4k;	//Для вычисления k4 из Рунге-Кутты
	k4x = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);
	k4k = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * dimArr);

	double x;
	for (i=0; i<dimArr; i++)		//Задание начальных условий
	{
		xCur = 2.0*M_PI*i*dx;
		x = i;
		//fx[i][0] = U0(alph*(x-dimArr/2), 10 , 0 ); 
		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]);
	}

	while (tCur < TIME)		// Старт цикла
	{
		tCur += t;
		//cout << tCur << endl;

		f_nonlinear(fk,NULL,dimArr,k1k);
		Aliasing(k1k);
		f_nonlinear(fk,k1k,dimArr,k2k);
		Aliasing(k2k);
		f_nonlinear(fk,k2k,dimArr,k3k);
		Aliasing(k3k);
		for (i = 0; i < dimArr; i++)
		{
			k3k[i][0] = k3k[i][0]*2;
			k3k[i][1] = k3k[i][1]*2;
		}
		f_nonlinear(fk,k3k,dimArr,k4k);
		Aliasing(k4k);


		//Пересчет РК 4 порядка
		for (i=0; i<dimArr; i++)		
		{
			xCur = 2.0*M_PI*i*dx;
			fk[i][0] = fk[i][0] + t*(k1k[i][0]+2*k2k[i][0] + k3k[i][0] + k4k[i][0])/6.0;
			fk[i][1] = fk[i][1] + t*(k1k[i][1]+2*k2k[i][1] + k3k[i][1] + k4k[i][1])/6.0;
		}
		Aliasing(fk);
	}

	//Переход в обычное пространство
	fftw_execute(fk_2_fx_); 
		
	for (i=0; i<dimArr; i++)		
	{
		xCur = 2.0*M_PI*i*dx;
		fprintf(Euler,"%1.5f   %1.5f   %1.5f\n",xCur, fx[i][0],fx[i][1]);
	}


	fftw_destroy_plan(fx_2_fk);
	fftw_destroy_plan(fk_2_fx_);
	fftw_free(fx);
	fftw_free(fk);



	
	fftw_free(k1x);
	fftw_free(k1k);

	fclose(ini);
	fclose(res);
	fclose(fk_data);
	fclose(Euler);


}