/work_files/curs/curs/parabolicVector_solver.cpp

  1#include "parabolicVector_solver.h"
  2
  3double **CParabolicVector_solver:: getTable(){
  4	return pTable;
  5}
  6
  7int CParabolicVector_solver::getK(){
  8	return K;
  9}
 10
 11int CParabolicVector_solver::getN(){
 12	return N;
 13}
 14
 15double * CParabolicVector_solver::getRightBorder()
 16{
 17	double *temp_arr = new double[K];
 18	for (int k =0; k < K; k++)
 19	{
 20		temp_arr[k] = pTable[k][N-1];
 21	}
 22	return temp_arr;
 23}
 24
 25CParabolicVector_solver::CParabolicVector_solver()
 26{
 27
 28}
 29
 30CParabolicVector_solver::~CParabolicVector_solver(){
 31	if (pTable != NULL)
 32	{
 33		for (int i = 0; i < K; ++i)
 34		{
 35			delete [] pTable[i];
 36		}
 37		delete [] pTable;
 38	pTable = NULL;
 39	}
 40	delete[] q_arr;
 41	q_arr = NULL;
 42
 43}
 44
 45void CParabolicVector_solver::setGrid(int _N, int _K, double _xb, double _xe, double _maxT)
 46{
 47	stepT = _maxT / _K;
 48	stepX = (_xe - _xb) / _N;
 49
 50	if (pTable != NULL)
 51	{
 52		for (int i = 0; i < K; ++i)
 53		{
 54			delete [] pTable[i];
 55		}
 56		delete [] pTable;
 57		pTable = NULL;
 58	}
 59
 60
 61	N = _N + 1;
 62	K = _K + 1;
 63	xb = _xb;
 64	xe = _xe;
 65	maxT = _maxT;
 66
 67//
 68
 69	pTable = new double* [K];
 70	q_arr = new double[K];
 71	for (int i = 0; i < K; ++i)
 72	{
 73		pTable[i] = new double[N];
 74		q_arr[i] = 0.0;
 75	}
 76	
 77}
 78
 79void CParabolicVector_solver::setTask(	double _a, double _b, double _c, double(*_f)(double, double),
 80							   double _alpha1, double _betta1, double _alpha2, double _betta2,
 81							   double * _left_border, double (*_phi2)(double, double), 
 82							   double (*_psi)(double, double)) // Only after setGrid;
 83{
 84	a = _a;
 85	b = _b;
 86	c = _c;
 87	f = _f;
 88	alpha1 = _alpha1;
 89	alpha2 = _alpha2;
 90	betta1 = _betta1;
 91	betta2 = _betta2;
 92	phi2 = _phi2;
 93	psi = _psi;
 94	setBorder(_left_border);
 95	
 96}
 97
 98void CParabolicVector_solver::solve_implicit()
 99{
100	double sigma = (a * a * stepT) / (stepX * stepX);
101
102	//????????? ?????? ?????? ?? ???????	psi(x, 0) = sin(x)
103	for (int i = 0; i < N; ++i)
104	{
105		pTable[0][i] = psi(xb + i * stepX, 0.0);
106	}
107
108	//??????? ??? ?????? ????????
109	double* A = new double[N];
110	double* B = new double[N];
111	double* C = new double[N];
112	double* D = new double[N];
113	double* P = new double[N];
114	double* Q = new double[N];
115	double* x = new double[N];
116
117	double csigma = a * a * stepT / (stepX * stepX); // only constant, don't want to copy-paste
118	double csigma1 = 2.0 * a * a * alpha1 / (2.0 * a * a * stepX - b * stepX * stepX);	// only constant, don't want to copy-paste
119	double csigma2 = 2.0 * a * a * alpha2 / (2.0 * a * a * stepX + b * stepX * stepX);	// only constant, don't want to copy-paste
120
121	for (int i = 1; i < K; ++i)
122	{
123		A[0] = 0.0;
124		B[0] = csigma1 * ( ((c * stepX * stepX) / (2.0 * a * a)) - 1.0 - ((stepX * stepX) / (2.0 * a * a * stepT))) + betta1;
125		C[0] = csigma1;
126		D[0] = phi0(xb, i ) - pTable[i-1][0] * (csigma1 * ((stepX * stepX) / (2.0 * a * a * stepT))) - csigma1 * ((stepX * stepX) / (2.0 * a * a)) * f(xb, i * stepT);
127	// ??? ? ?????? ???? - ??????????? i -?? ??????? ?? ???????.
128		for (int j = 1; j < N - 1; ++j)
129		{
130			A[j] = csigma - ((b * stepT) / (2.0 * stepX));
131			B[j] = (-2.0 * csigma) - 1.0 + c * stepT;
132			C[j] = csigma + ((b * stepT) / (2.0 * stepX));
133			D[j] = -pTable[i-1][j] - f(xb + j * stepX, i * stepT) * stepT;
134
135			if (fabs(B[i]) < fabs(A[i]) + fabs(C[i]))
136			{
137				printf("error in coeff\n");
138			}
139		}
140
141		A[N - 1] = -csigma2;
142		B[N - 1] = csigma2 * (1.0 + ((stepX * stepX) / (2.0 * a * a * stepT)) - ((c * stepX * stepX) / (2.0 * a * a))) + betta2;
143		C[N - 1] = 0.0;
144		D[N - 1] = phi2(xe, i * stepT) + pTable[i-1][N - 1] * (csigma2 * stepX * stepX / (2.0 * a * a * stepT)) + f(xe, i * stepT) * csigma2 * stepX * stepX / (2.0 * a * a);
145
146		// ?????? ????????
147		double aa, bb, cc;
148
149		bb = B[0];
150		cc = C[0];
151
152		P[0] = -cc / bb;
153		Q[0] = D[0] / bb;
154
155		for (int k = 1; k < N - 1; ++k)
156		{
157			aa = A[k];
158			bb = B[k];
159			cc = C[k];
160			P[k] = -cc / (bb + aa*P[k-1]);
161			Q[k] = (D[k] - aa*Q[k-1]) / (bb + aa*P[k-1]);
162		}
163
164		aa = A[N - 1];
165		bb = B[N - 1];
166		P[N - 1] = 0;
167		Q[N - 1] = (D[N - 1] - aa * Q[N - 2]) / (bb + aa * P[N - 2]);
168
169		x[N - 1] = Q[N - 1];
170		for (int k = N - 2; k >= 0; k--)
171		{
172			x[k] = P[k] * x[k+1] + Q[k];
173		}
174		//????? ????????
175
176		for (int j = 0; j < N; ++j)
177			pTable[i][j] = x[j];
178	}
179
180	delete [] A;
181	delete [] B;
182	delete [] C;
183	delete [] D;
184	delete [] P;
185	delete [] Q;
186	delete [] x;
187}
188
189double CParabolicVector_solver::phi0(double x, int i){
190	return -q_arr[i];
191}
192
193void CParabolicVector_solver::setBorder(double *_arr){
194	for (int i = 0; i < K; i++){
195		q_arr[i] = _arr[i];
196	}
197}
198