/project/curs/curs/parabolicVector_solver.cpp

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