test_fieldeq_sc.cpp

#include <iostream>
#include <cmath>
#include <stdio.h>
using namespace std;

double curtime;
int size;
double eps = 1.0;
double k = 0.001;
double dt = 0.1;
double cspeed = 1.0;
double r = cspeed*dt;
double omega = cspeed*k/sqrt(eps);

void set_initial(double * E_prev, double * E_prev_prev)
{
	for (int i = 0; i < size; i++){
		E_prev_prev[i] = 1.0;
		E_prev[i] = 1.0;
	}
	return;
}

void massive_set(double * a, double * b, double * c, double * f,double * E_prev, double * E_prev_prev)		//Неявная схема
{
	
	for(int i = 1; i < size-1; i++)
		{
			*(a+i) = sin(-k*i-omega*curtime)-4*k*cos(-k*i-omega*curtime); //4*k*cos(-k*i-omega*curtime)
			*(b+i) = sin(-k*i-omega*curtime)*2.0*(1.0+2.0*(2*eps-eps)/(r*r)) - 8.0*omega*eps/(r*r)*dt*cos(-k*i-omega*curtime);// - 8.0*k*cos(-k*i-omega*curtime);
			*(c+i) = sin(-k*i-omega*curtime)+4*k*cos(-k*i-omega*curtime);
			*(f+i) = sin(-k*i-omega*curtime)*2.0*(E_prev[i+1] - 2.0*E_prev[i] + E_prev[i-1] + 4.0*eps/(r*r)*E_prev[i])
					+ sin(-k*i-omega*curtime)*(E_prev_prev[i+1] - 2.0*E_prev_prev[i] + E_prev_prev[i-1] - 4.0*eps/(r*r)*E_prev_prev[i])
					- cos(-k*i-omega*curtime)*4.0*eps*dt*2.0*omega/(r*r)*E_prev[i];
		}
	
	
	//Зададим условие на значение произодной слева (T'=Ua, x=0) и значение функции справа (T=Ub, x=1)


	//Set left border conditions: (derivat.) - задание условия на производную слева
	*a = 0.0;
	*b = 1.0;
	*c = 1.0;
	*f = 0.0;
	
	/*
	// Задание условия на значение слева, а не на производную
	*b = 1.0;
	*c = 0.0;
	*f = Ua;
	*/

	//Задание условия на значения справа
	*(a+size-1) = 0.0;
	*(b+size-1) = 1.0;
	*(c+size-1) = 0.0;
	*(f+size-1) = 1.0;
		
	return;
}

void massive_get(double * a, double * b, double * c, double * f, double * beta, double * z)		//Заполнение массивов бета, альфа и Z
{
	*beta =  *c / *b;
	*z = *f / *b;

	for(int i = 1; i < size; i++)
	{
		//if (abs(*(b+i) - (*(a+i))*(*(beta+i-1))) < 0.001)
		//	cout << (*(b+i) - (*(a+i))*(*(beta+i-1))) << endl;
		
		*(beta+i) = *(c+i) / (*(b+i) - (*(a+i))*(*(beta+i-1)));

		//if (i == size-1)
		//	cout << "a: " << *(a+i) << "b: " << *(b+i) << "f: " << *(f+i) << "z-1: " << *(z+i-1) << "beta-1: " << *(beta+i-1) << endl;

		*(z+i) = (*(f+i) +(*(a+i))*(*(z+i-1)))/(*(b+i) - (*(a+i))*(*(beta+i-1)));
	}

	//*(beta+size-1) = 0.0;
	//*(z+size-1) = 1.0;
	return;
}

void get_solution( double * beta, double * z,double * E_new)		//конечный пересчет
{
	*(E_new+size-1) = 1.0;

	for (int i = size-2; i >= 0; i--){
		*(E_new+i) = *(beta+i) * (*(E_new+i+1)) + *(z+i);
		
		//if (*(E_new+i) > 2.0)
			//cout << *(E_new+i+1) << "  " <<  i << "  " <<  curtime << endl;

		//if (curtime > 295) cout <<  *(beta+i) << " :beta, " << *(z+i) << " :z;" << " curtime: " << curtime << " position: " << i << endl;
	}
		
	return;
}

int main()
{
	FILE * file = fopen("file.txt","w");
		FILE * filea = fopen("filea.txt","w");
		FILE * fileb = fopen("fileb.txt","w");
		FILE * filec = fopen("filec.txt","w");
		FILE * filef = fopen("filef.txt","w");
		FILE * filebeta = fopen("filebeta.txt","w");
		FILE * filez = fopen("filez.txt","w");

	size = 6280;
	curtime = 0;	
	double * E_prev_prev = new double [size * sizeof(double)];
	double * E_prev = new double [size * sizeof(double)];
	double * E_new = new double [size * sizeof(double)];

	double * a = new double[size * sizeof(double)];	//Подсчет коэффициентов a,b,c из схемы и beta,z -  
	double * b = new double[size * sizeof(double)]; //- промежуточных коэффициентов в неявной схеме
	double * c = new double[size * sizeof(double)];
	double * f = new double[size * sizeof(double)];
 
	double * beta = new double[size * sizeof(double)];
	double * z = new double[size * sizeof(double)];

	double final_time = 1100.0;

	set_initial(E_prev,E_prev_prev);

	//for (int i = 0; i < size; i++)
	//	fprintf(file, "%e\n", E_prev_prev[i]);

	while (curtime < final_time)
	{
		curtime = curtime + dt;

		massive_set(a,b,c,f,E_prev,E_prev_prev);

		massive_get(a,b,c,f,beta,z);

		get_solution(beta,z,E_new);

		E_prev_prev = E_prev;

		E_prev = E_new;

		

		/*for (int i = 0; i < size; i++){
			fprintf(file, "%e\n", *(E_new+i)*sin(-k*i-omega*curtime));
			fprintf(filea, "%e\n", *(a+i));
			fprintf(fileb, "%e\n", *(b+i));
			fprintf(filec, "%e\n", *(c+i));
			fprintf(filef, "%e\n", *(f+i));
			fprintf(filebeta, "%e\n", *(beta+i));
			fprintf(filez, "%e\n", *(z+i));
		}*/
		
		

		//cout << curtime << endl;
		//cin.get();

		//fprintf(file, "%e\n", *(E_new+400));
		}

		for (int i = 0; i < size; i++){
			fprintf(file, "%e\n", *(E_new+i)*sin(-k*i-omega*curtime));
			fprintf(filea, "%e\n", *(a+i));
			fprintf(fileb, "%e\n", *(b+i));
			fprintf(filec, "%e\n", *(c+i));
			fprintf(filef, "%e\n", *(f+i));
			fprintf(filebeta, "%e\n", *(beta+i));
			fprintf(filez, "%e\n", *(z+i));
		}

		fclose(file);
		fclose(filea);
		fclose(fileb);
		fclose(filec);
		fclose(filef);
		fclose(filebeta);
		fclose(filez);



	return 0;
}