#include "stdafx.h"
#define _USE_MATH_DEFINES
# include <iostream>
#include <fstream>
#include <cmath>
#include <math.h>
using namespace std;

const double D = 1.0;
const double lambda = 1; //dla metod posrednich
double h = 0.02; //krok dla x
double dt = (lambda*h*h) / D; //krok dla t
const double x_max = 1; //maksymalny przedzial x
const double x_min = 0; //minimalny przedzial x
const double t_max = 0.5; //maksymalny przedzial x
const double t_min = 0; //minimalny przedzial x

//wymiary macierzy do obliczeń numerycznych:
int M = (int)((x_max - x_min) / h) + 1; //liczba kolumn (x) 11
int N = (int)((t_max / dt) + 1); //liczba wierszy(t) 51

void rozwiazanie_analityczne(double **Analityczne)
{
	fstream zapis;
	double x = x_min + h;	//zmienna wartosci przestrzennej
	double t = dt;			//zmienna kroku czasowego

	zapis.open("rozw_analit.txt", ios::out);

	//wyzerowanie macierzy rozwiązań analitycznych
	for (int i = 0; i < N; i++)
		for (int j = 0; j < M; j++)
			Analityczne[i][j] = 0.0;
	//wpisanie do macierzy wartości wynikłych z warunków początkowych i brzegowych
	//warunek poczatkowy
	for (int j = 0; j < M; j++)
	{
		x = h*(double)j + x_min;
		Analityczne[0][j] = 1.0 + cos(M_PI*x);
		x = x + h;
	}
	//warunki brzegowe
	for (int i = 1; i < N; i++)
	{
		Analityczne[i][0] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI * 0);
		Analityczne[i][M - 1] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI * 1); // zapisujemy do macierzy wartości wyliczone
		t = t + dt;											//krok t
	}
	t = dt;
	x = x_min + h;
	//POZOSTAŁE PUNKTY SIATKI
	for (int i = 1; i < N; i++)
	{
		for (int j = 1; j < M - 1; j++)
		{
			Analityczne[i][j] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI*x);	// zapisujemy do macierzy wartości wyliczone
			x = x + h;
		}
		x = x_min + h;	//przywrócenie x do początku
		t = t + dt;
	}

	//zapisujemy wyniki rozwiązania analitycznego do pliku
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < M; j++)
		{
			zapis << Analityczne[i][j] << "; ";
		}
		zapis << endl;
	}
	zapis.close();
}

void operacje_na_macierzy(double *l, double *u, double *d)
{
	for (int i = 1; i < M; i++)
	{
		d[i] = d[i] - l[i - 1] * u[i - 1] / d[i - 1];  //faza eliminacji wprzód
	}
}

void operacje_na_wektorze(double *b, double *u, double *d, double *l, double *x)
{
	for (int i = 1; i < M; i++)
	{
		b[i] = b[i] - l[i - 1] * b[i - 1] / d[i - 1];
	}
	x[M - 1] = b[M - 1] / d[M - 1];
	for (int i = M - 2; i >= 0; i--)
	{
		x[i] = (b[i] - u[i] * x[i + 1]) / d[i];
	}
}

void SOR(double **A, double *b, double *x)
{
	double *x0_pom = new double[M];
	double *x0 = new double[M];
	double sum = 0.0;
	double EST = 1e-6;
	double omega = 0.4;
	double max_diff;
	double diff;
	for (int i = 0; i < M; i++)
		x0[i] = 0.0;

	for (int i = 0; i < 100; i++)  //liczba iteracji
	{
		for (int x = 0; x < M; x++)
			x0_pom[x] = x0[x];

		for (int j = 0; j < M; j++)
		{
			sum = (1.0 - 1.0 / omega)*A[j][j] * x0[j];
			for (int k = j + 1; k <M; k++)
				sum += A[j][k] * x0_pom[k];

			for (int k = 0; k < j; k++)
				sum += A[j][k] * x0_pom[k];

			x0[j] = (b[j] - sum)*omega / A[j][j];
		}
		//sprawdzenie czy różnica pomiędzy kolejnymi przybliżeniami jest mniejsza od zadanej wielkości
		//poszukiwanie normy max z wektora rozwiązań
		max_diff = fabs(x0_pom[0] - x0[0]);
		for (int k = 1; k < M; k++)
		{
			diff = fabs(x0_pom[k] - x0[k]);
			if (diff > max_diff)
				max_diff = diff;
		}
		if (max_diff<EST)
			break;
	}
	for (int i = 0; i < M; i++)
		x[i] = x0[i];
	delete[] x0_pom;
	delete[] x0;
}

void Crank_Nicolson(double **A)
{
	fstream zapis;
	double *l, *d, *u;			// wektory L, D, U
	double *b, *xrozw;			// wektor wyrazów wolnych oraz wektor rozwiązań
	double xpom;
	l = new double[M];
	d = new double[M];
	u = new double[M];
	b = new double[M];
	xrozw = new double[M];

	zapis.open("thomas.txt", ios::out);
	//wpisanie do macierzy wartości wynikłych z warunków początkowych i brzegowych
	//WARUNEK POCZĄTKOWY
	for (int j = 0; j < M; j++)
	{
		xpom = h*(double)j + x_min;
		A[0][j] = 1.0 + cos(M_PI*xpom);
	}

	double t = dt;
	//BRZEGOWE
	for (int i = 1; i < N; i++)
	{
			A[i][0] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI*0);
			A[i][M -1] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI*1); // zapisujemy do macierzy wartości wyliczone
			t = t + dt;											//krok t
	}
	t = dt;
	// wypełnienie l, d, u
	for (int j = 0; j < M; j++)
		l[j] = lambda / 2.0;		//LOWER
	l[M - 2] = -1/h;   //M-2

	for (int j = 1; j < M - 1; j++)		//DIAGONAL
		d[j] = -(1.0 + lambda);
	d[0] = -1.0 / h;
	d[M - 1] = 1.0/h;

	for (int j = 1; j < M; j++)	//UPPER
		u[j] = lambda / 2.0;
	u[0] = 1.0/h;

	//operacje na l,d,u
	operacje_na_macierzy(l, u, d);

	//pętla do obliczania przybliżeń dla kroku czasowego
	for (int i = 0; i < N - 1; i++)
	{
		//wypełniamy wektor wyrazów wolnych
		b[0] = 0.0;
		b[M - 1] = 0.0;

		for (int j = 1; j < M - 1; j++)
		{
			b[j] = -(lambda / 2.0)*A[i][j - 1]  - (1 - lambda)*A[i][j]  -(lambda / 2.0)*A[i][j + 1];
		}


		//druga część Thomasa
		operacje_na_wektorze(b, u, d, l, xrozw);

		//przepisanie wyników do macierzy A
		for (int j = 1; j<M - 1; j++)
			A[i + 1][j] = xrozw[j];
	}

	//zapis do pliku macierzy A
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < M; j++)
		{
			zapis << A[i][j] << "; ";
		}
		zapis << endl;
	}
	zapis.close();
	delete[] l;
	delete[] d;
	delete[] u;
	delete[] b;
	delete[] xrozw;
}

void Crank_Nicolson_SOR(double **A)
{
	fstream zapis;
	double ** macierz; //MACIERZ POMOCNICZA
	double *b, *xrozw; //WEKTOR WYRAZÓW WOLNYCH ORAZ POMOCNICZY WKTOR ROZWIĄZAŃ
	double xpom;

	b = new double[M];
	xrozw = new double[M];
	macierz = new double *[M];
	for (int i = 0; i < M; i++)
		macierz[i] = new double[M];
	zapis.open("SOR.txt", ios::out);
	//wpisanie do macierzy wartości wynikłych z warunków początkowych i brzegowych
	//WARUNEK POCZĄTKOWY

	for (int j = 0; j < M; j++)
	{
		xpom = h*(double)j + x_min;
		A[0][j] = 1.0 + cos(M_PI*xpom);
		xpom = xpom + h;
	}

	double t = dt;
	//BRZEGOWE
	for (int i = 1; i < N; i++)
	{
		A[i][0] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI * 0);
		A[i][M - 1] = 1 + exp(-(M_PI*M_PI)*D*t)*cos(M_PI * 1); // zapisujemy do macierzy wartości wyliczone
		t = t + dt;											//krok t
	}
	t = dt;

	//WYZEROWANIE MACIERZY POMOCNICZEJ
	for (int i = 0; i < M; i++)
		for (int j = 0; j < M; j++)
			macierz[i][j] = 0.0;
	// PRZYGOTOWANIE MACIERZY TRÓJDIAGONALNEJ
	macierz[1][0] = lambda / 2.0;
	macierz[0][0] = -1.0/h;
	macierz[0][1] = 1.0/h;
	for (int j = 1; j < M - 1; j++)
	{
		macierz[j + 1][j] = lambda / 2.0;		//LOWER
		macierz[j][j] = -(1.0 + lambda);		//DIAGONAL
		macierz[j][j + 1] = lambda / 2.0;		//UPPER
	}
	macierz[M - 1][M - 1] = 1.0/h;
	macierz[M - 1][M - 2] = -1.0/h;


	//ROZPOCZYNAMY PĘTLE W KTÓREJ OBLICZAMY KOLEJNE PRZYBLIŻENIA DLA KOLEJNYCH
	//PUNKTÓW NA LINII CZASOWEJ
	for (int i = 0; i < N - 1; i++)
	{

		//WYPEŁNIAMY WARTOŚCIAMI WEKTOR WYRAZÓW WOLNYCH
		b[0] = 0.0;
		b[M - 1] = 0.0;
		for (int j = 1; j < M - 1; j++)
			b[j] = -(lambda / 2.0)*A[i][j - 1] - (1 - lambda)*A[i][j] - (lambda / 2.0)*A[i][j + 1];


		//ROZWIĄZUJEMY UKŁAD ALGEBRAICZNYCH RÓWNAŃ LINIOWYCH
		//METODĄ SOR
		SOR(macierz, b, xrozw);

		//PRZEPISUJEMY WARTOŚCI Z WEKTORA DO MACIERZY ROZWIĄZAŃ
		for (int j = 1; j<M - 1; j++)
			A[i + 1][j] = xrozw[j];
	}
	//ZAPISUJEMY MACIERZ ROZWIĄZAŃ NUMERYCZNYCH METODĄ SOR DO PLIKU
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < M; j++)
		{
			zapis << A[i][j] << "; ";
		}
		zapis << endl;
	}
	zapis.close();

	for (int i = 0; i < M; i++)
		delete[] macierz[i];
	delete[] macierz;
	delete[] b;
	delete[] xrozw;
}

double blad(double **Blad, double **rozwiazanie, double **analityczne, int met)
{
	double *blad_max = new double[N];

	fstream zapis;
	fstream zapis2;
	if (met)	//SPRAWDZAMY DLA KTÓREJ METODY OBLICZAMY BŁĄD
	{//I OTWIERAMY ODPOWIEDNIE PLIKI DO ZAPISU
		zapis.open("blad_max_thomas.txt", ios::out);
		zapis2.open("blad_thomas.txt", ios::out);
	}
	else
	{
		zapis.open("blad_max_sor.txt", ios::out);
		zapis2.open("blad_SOR.txt", ios::out);
	}
	// ZERUJEMY WEKTOR BŁĘDÓW MAX
	for (int i = 0; i < N; i++)
		blad_max[i] = 0.0;
	//OBLICZAMY MACIERZ BŁĘDÓW DLA ZADANEJ METODY
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < M; j++)
		{
			Blad[i][j] = fabs(analityczne[i][j] - rozwiazanie[i][j]);
			zapis2 << Blad[i][j] << "; ";
			//OBLICZAMY BŁĄD MAX DLA DANEJ CHWILI CZASU
			if (Blad[i][j]>blad_max[i])
				blad_max[i] = Blad[i][j];
		}
		zapis2 << endl;
	}
	for (int i = 0; i < N; i++)
		zapis << blad_max[i] << " ";

	//OBLICZAMY MAX WARTOŚĆ BEZWZGLĘDNĄ DLA CZASU TMAX
	double tmp = Blad[N - 1][0];
	for (int i = 1; i < M; i++)
	{
		if (Blad[N - 1][i]>tmp)
			tmp = Blad[N - 1][i];
	}
	if (met)
		cout << "Thomas:" << h << " " << tmp <<endl;
	else
		cout << "SOR " << h << " " << tmp << endl;
	return tmp;
	zapis.close();
	zapis2.close();
}

void zapisz_kroki()
{
	double czas = 0.0;
	double przestrzen;
	fstream k_czas;
	fstream k_przestrzenne;

	k_czas.open("kroki_czasowe.txt", ios::out);
	k_przestrzenne.open("kroki_przestrzenne.txt", ios::out);

	przestrzen = x_min;
	//ZAPISUJEMY DO PLIKU KOLEJNE PUNKTY CZASOWE
	for (int i = 0; i < N; i++)
	{
		czas = dt*(double)i;
		k_czas << czas << endl;
	}
	//ZAPISUJEMY DO PLIKU KOLEJNE PUNKTY PRZESTRZENNE
	for (int i = 0; i < M; i++)
	{
		przestrzen = h*(double)i + x_min;
		k_przestrzenne << przestrzen << " ";
	}
	k_czas.close();
	k_przestrzenne.close();
}

int _tmain(int argc, _TCHAR* argv[])
{
	fstream bladd;
	bladd.open("bladmaxkrok.txt", ios::out);
	int tmp = 6;
	while (tmp<7)
	{
		dt = (lambda*h*h) / D;
		N = (int)(t_max / dt) + 2;
		M = (int)((x_max - x_min) / h) + 1;
		double **Analityczne;			//macierz rozwiązań analitycznych
		double **rozw_thomas;//macierz rozwiązań numerycznych algorytmem thomasa
		double **rozw_SOR;		//macierz rozwiązań numerycznych metodą SOR
		double **Blad;
		rozw_thomas = new double*[N];
		rozw_SOR = new double*[N];
		Analityczne = new double *[N];
		Blad = new double *[N];

		for (int i = 0; i < N; i++)
		{											//ALOKACJA PAMIĘCI
			Analityczne[i] = new double[M];
			rozw_thomas[i] = new double[M];
			rozw_SOR[i] = new double[M];
			Blad[i] = new double[M];
		}

		bladd << h << " ";
		zapisz_kroki();										//ZAPISUJEMY KROKI CZASOWE I PRZESTRZENNE DO PLIKU
		rozwiazanie_analityczne(Analityczne);				//OBLICZAMY I ZAPISUJEMY ROZWIĄZANIA ANALITYCZNE
		Crank_Nicolson(rozw_thomas);						//OBLICZAMY ROZWIĄZANIA NUMERYCZNE ZA POMOCĄ ALGORYTMU THOMASA
		bladd << blad(Blad, rozw_thomas, Analityczne, 1) << " " << endl;//OBLCZAMY BŁĘDY ORAZ ZAPISUJEMY JE DO PLIKU
		Crank_Nicolson_SOR(rozw_SOR);						//OBLICZAMY ROZWIĄZANIA NUMERYCZNE ZA POMOCĄ METODY SOR
		bladd << blad(Blad, rozw_SOR, Analityczne, 0) << endl;	//OBLCZAMY BŁĘDY ORAZ ZAPISUJEMY JE DO PLIKU

		for (int i = 0; i < N; i++)		//ZWALNIAMY PAMIĘĆ
		{
			delete[] Analityczne[i];
			delete[] rozw_thomas[i];
			delete[] rozw_SOR[i];
			delete[] Blad[i];
		}
		delete[] Analityczne;
		delete[] rozw_thomas;
		delete[] rozw_SOR;
		delete[] Blad;
		h = h/2;
		tmp++;
	}
	bladd.close();
	system("pause");
	return 0;
}