Теплопроводность неявная + Кранк-николсон

gistfile1.cpp

#define _CRT_SECURE_NO_DEPRECATE
#include "stdio.h"
#include "math.h"
#include <iostream>
using namespace std;

/*
Общие вопросы по теплопроводности:
Устойчивость для явной схемы и показать развал
Декремент затухания и сравнить с аналитикой

Общее решение T = Summ[ exp(-pi^2 n^2 t) * Bn * sin(pi*n*x)]; Bn задаются начальными условиями
*/

double lambda = 1.0; //Коэффициент теплопроводности

double t = 0.01; //Условие монотонности для схемы Кранка-Николсона t < h^2
double h = 0.2;

double a = 0.0;
double b = 1.0;

double Ua = 1.0; //left border condition
double Ub = 0.0; //right border condition

int size = 1+(b-a)/h;


double function(double x)	//T(t = 0)
{
	return x*x*(1-x);
	
}

void set_conditions(double * u)		//задание начальные условия
{
	for (int i = 0; i*h <= (b-a); i++)
		*(u+i) = function(i*h+a);
}

void massive_set(double * a, double * b, double * c, double * f,double * u)		//Неявная схема
{
	for(int i = 0; i < size-1; i++)
		{
			*(a+i) = -1.0/(h*h);
			*(b+i) = -2.0/(h*h)-1.0/t;
			*(c+i) = -1.0/(h*h);
			*(f+i) = -*(u+i)/t;
		}
	*a = 0.0;
	*(c+size-1) = 0.0;
	
	
	//Set left border conditions: (derivat.) - задание условия на производную слева
	*b = 2.0/3.0;
	*c = 2.0/3.0 - h*h/(3.0*t);
	*f = *(u+1)*h*h/(3.0*t) - 2.0/3.0 * h * Ua;
	
	/*
	//Set left border conditions:  - задание условия на значение слева
	*b = 1.0;
	*c = 0.0;
	*f = Ua;
	*/

	//Set right border condition
	*(a+size-1) = 0;
	*(b+size-1) = 1;
	*(f+size-1) = Ub;
		
	return;
}

void massive_set_KN(double * a, double * b, double * c, double * f,double * u)		//Схема Кранка-Николсона
{
	for(int i = 0; i < size-1; i++)
		{
			*(a+i) = -1.0/(2*h*h);
			*(b+i) = -1.0/(h*h)-1.0/t;
			*(c+i) = -1.0/(2*h*h);
			*(f+i) = -*(u+i)*(1/t-1/(h*h))-*(u+i+1)/(2*h*h)-*(u+i-1)/(2*h*h);
		}
	*a = 0.0;
	*(c+size-1) = 0.0;
	
	/*
	//Set left border conditions: (derivat.)
	*b = 2.0/3.0;
	*c = 2.0/3.0 - h*h/(3.0*t);
	*f = *(u+1)*h*h/(3.0*t) - 2.0/3.0 * h * Ua;
	*/
	//Set left border conditions:
	*b = 1.0;
	*c = 0.0;
	*f = Ua;

	//Set right border condition
	*(a+size-1) = 0;
	*(b+size-1) = 1;
	*(f+size-1) = Ub;
		
	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-1; i++)
	{
		*(beta+i) = *(c+i) / (*(b+i) - (*(a+i))*(*(beta+i-1)));
		*(z+i) = (*(f+i) +(*(a+i))*(*(z+i-1)))/(*(b+i) - (*(a+i))*(*(beta+i-1)));
	}
	return;
}

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

	for (int i = size-2; i >= 0; i--)
		*(u+i) = *(beta+i) * (*(u+i+1)) + *(z+i);

	return;
}

void main()
{
	FILE * file = fopen("file.txt", "w");
	
	double * u = new double[size * sizeof(double)];
	double * u_next = new double[size * sizeof(double)];
	
	set_conditions(u); //Set start temperature distribution
	
	//Massives for non-clear scheme
	double * a = new double[size * sizeof(double)];
	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 time = 0.1;
	double time_current = 0.0;
	
	while (time_current < time)
	{
		time_current = time_current + t;
		massive_set(a,b,c,f,u);
		massive_get(a,b,c,f,beta,z);
		get_solution(beta,z,u_next);
		for (int i = 0; i < size; i++)
			*(u+i) = *(u_next + i);
	}

	for (int i = 0; i < size; i++)
			fprintf(file, "%e\n",  *(u+i));//fprintf(file, "%e %e\n", i*h, *(u+i)); - в комменте вывод с координатами

	return;
}