PageRenderTime 64ms CodeModel.GetById 19ms app.highlight 41ms RepoModel.GetById 1ms app.codeStats 0ms

/project/curs/curs/backTime_solver.cpp

http://cm-kp-406.googlecode.com/
C++ | 284 lines | 167 code | 58 blank | 59 comment | 13 complexity | 97d420858fb20ed1d752c375ba1ee0ba MD5 | raw file
  1#include "backTime_solver.h"
  2
  3
  4CBackTimeSolve :: CBackTimeSolve()
  5{
  6	pTable = NULL;
  7	analytic = NULL;
  8	right_border = NULL;
  9}
 10
 11CBackTimeSolve :: ~CBackTimeSolve()
 12{
 13	clearArray();
 14	analytic = NULL;
 15	right_border = NULL;
 16}
 17
 18
 19double ** CBackTimeSolve::getTable()
 20{
 21	return pTable;
 22}
 23
 24
 25double * CBackTimeSolve::getGradient(double * prev_grad)
 26{
 27	for (int k =0; k < K; k++)
 28	{
 29		prev_grad[k] = -pTable[k][0];
 30	}
 31	return prev_grad;
 32
 33}
 34
 35
 36/*
 37void ReverseSolve::CheckSolve(string f_name)
 38{
 39	ofstream out(f_name.c_str());
 40	out << fixed;
 41	float temp;
 42	for (int k = 0; k < K; k++)
 43	{
 44		for (int i = 0; i < N; i++)
 45		{
 46			temp = r(stX + i*h, k*tau);
 47			out << setprecision (5) << "t = " << k*tau << "   x = " << stX + i*h << "   u(x,y) = " << u_array[k][i] << "   u*(x,y) = " << temp << "  Difference = " << fabs(u_array[k][i] - temp) << endl;
 48		}
 49		out << endl << endl;
 50	}
 51}*/
 52
 53
 54
 55
 56
 57void CBackTimeSolve:: clearArray()
 58{
 59	if (pTable != NULL)
 60	{
 61		for (int k = 0; k < K; k++)
 62		{
 63			delete [] pTable[k];
 64		}
 65		delete [] pTable;
 66	}
 67	pTable = NULL;
 68}
 69
 70
 71/*
 72void ReverseSolve::CheckResult(double **correct, string f_name)
 73{
 74	ofstream out(f_name.c_str());
 75	out << fixed;
 76	double temp;
 77	for (int k = 0; k < K; k++)
 78	{
 79		for (int i = 0; i < N; i++)
 80		{
 81			temp = fabs(correct[k][i] - u_array[k][i]);
 82			out << setprecision (5) << "t = " << k*tau << "   x = " << stX + i*h << "   u(x,y) = " << u_array[k][i] << "   u*(x,y) = " << correct[k][i] << "  Difference = " << temp << endl;
 83		}
 84		out << endl << endl;
 85	}
 86}*/
 87
 88
 89
 90
 91
 92
 93void CBackTimeSolve::checkResult( string f_name)
 94{
 95	ofstream out(f_name.c_str());
 96	out << fixed;
 97	for (int k = K-1; k >= 0; k--)
 98	{
 99		for (int i = 0; i < N; i++)
100		{
101			out << setprecision (5) << "t = " << T - k*tau << "   x = " << stX + i*h << "   u(x,y) = " << pTable[k][i] << endl;
102		}
103		out << endl << endl;
104	}
105}
106
107void CBackTimeSolve::setTask(double _a, double _b, double _c, 
108						   double _alpha, double _beta, double _gamma, double _delta, 
109						    double _stX, double _endX, double _Tr, int _N, int _K, 
110							double *_right_border, double *_analitic)
111{
112	clearArray();
113	a = _a;
114	b = _b;
115	c = _c;
116	alpha = _alpha;
117	beta = _beta;
118	gamma = _gamma;
119	delta = _delta;
120	l = _endX;
121	T = _Tr;
122	stX = _stX;
123	N = _N + 1;
124	K = _K + 1;
125	h = _endX/_N;
126	tau = _Tr/_K;
127	pTable = createTable();
128	right_border = _right_border;
129	analytic = _analitic;
130}
131
132double CBackTimeSolve::r(double x, double t)
133{
134	    return exp(-1*a*t)*sin(x);
135	//return sin(t) * cos(x);
136	//return exp(- a * a * t) * cos(x + b * t);
137	//return exp((c - a * a) * t) * sin(x + b * t);
138}
139
140double ** CBackTimeSolve::createTable()
141{
142	double ** new_array = new double*[K];
143	for (int i =0; i < K ; i++)
144	{
145		new_array[i] = new double[N];
146	}
147	return new_array;
148}
149
150void CBackTimeSolve::UseSweepMethod(double *d, double *A, double *B, double *C, int N, double *res)
151{
152	double* P = new double[N];
153	double* Q = new double[N];
154	double a,b,c;	
155	b = B[0];
156	c = C[0];
157	P[0] = -c / b;
158	Q[0] = d[0] / b;
159
160	for (int i=1; i < N-1; i++)
161	{
162		a = A[i];
163		b = B[i];
164		c = C[i];
165		P[i] = -c / (b + a*P[i-1]);
166		Q[i] = (d[i] - a*Q[i-1]) / (b + a*P[i-1]);
167	}
168
169	a = A[N-1];
170	b = B[N-1];
171	P[N-1] = 0;
172	Q[N-1] = (d[N-1] - a * Q[N-2]) / (b + a * P[N-2]);
173
174	res[N-1] = Q[N-1];					//SweapMethod
175	for (int i=N-2; i >= 0; i--)
176	{
177		res[i] = P[i]*res[i+1] + Q[i];
178	}
179}
180
181
182
183double CBackTimeSolve:: psiT(double x, double t)
184{
185	return 0.0;
186}
187
188// ????????????
189//double CBackTimeSolve:: psiT(double x, double t)
190//{
191//	return sin (x);
192//}
193
194
195
196double CBackTimeSolve:: f(double x, double t)
197{
198	return 0.0;
199}
200
201double CBackTimeSolve:: phi0(double x, double t)
202{
203	return 0.0;
204}
205
206
207double CBackTimeSolve:: phiL(double x, int k)
208{ 
209	return (2.0 * ( right_border[k] - analytic[k])) ; 
210	// ?? ??????? ?????? ?????? ??? ????????? ?? ?????? ??????? ??????? 
211	// ????????????? ??????? ????????(??????????????) ??????  
212}
213
214// ????????????
215//double CBackTimeSolve:: phiL(double x, double t)
216//{
217//	//return (2.0 * ( border[k] - tStar[k])); 
218//	return -1*exp(-1*a*t);
219//	//return -sin(t);
220//	//return (  exp(-a * a * t) * (cos(b * t) + sin (b*t))  );
221//	//return (  -exp((c - a*a)*t) * (cos (b * t) + sin (b * t))  );
222//}
223
224
225
226/*
227void CBackTimeSolve::setBorder(double * _right_border)
228{
229	right_border = _right_border;
230}
231*/
232
233void CBackTimeSolve::getSolve()
234{
235
236
237	double * A = new double[N];
238	double * B = new double[N];
239	double * C = new double[N];
240	double * d = new double[N];
241
242	double temp_a, temp_b, temp_c;
243	temp_a = (a/h/h - b/h/2.0)*tau;
244	temp_b = -(2*a*tau/h/h + 1 - c*tau);
245	temp_c = (a/h/h + b/2/h)*tau;
246
247	A[0] = 0;
248	B[0]= ((-3 + temp_a/temp_c )*alpha + 2* beta*h);
249	C[0]= (4 + temp_b/temp_c)*alpha;
250
251	A[N-1] = - (4 + temp_b/temp_a)*gamma;
252	B[N-1] = (3 - temp_c/temp_a)*gamma + 2*h*delta;
253	C[N-1] = 0;
254
255	for (int i = 1; i < N-1; i++)
256	{
257		A[i] = temp_a;
258		B[i] = temp_b;
259		C[i] = temp_c;
260	}
261
262	for (int i = 0; i < N; i++)
263	{
264		pTable[K-1][i] = psiT(stX + i*h, 0 ); //psiT(i ,0);
265	}
266
267	for (int k = K-2; k >= 0; k-- )
268	{
269		for (int i = 1; i < N; i++ )
270		{
271
272			d[i] = -pTable[k+1][i] - f(stX + i* h, ( k*tau) )*tau;
273
274		}
275
276		d[0] = phi0(0, ( k*tau) )*2*h + ( d[1] + f(0, (k*tau)) )*alpha/temp_c ;
277		d[N-1] = phiL(l, k )*2*h - (d[N-2] + f(l, k*tau) )*gamma/temp_a;
278		// ??? ? ?????? ???? ??????? k- ?? ???????? ?? ???????.
279		
280		//d[N-1] = phiL(l,( k* tau) )*2*h - (d[N-2] + f(l,( k*tau)) )*gamma/temp_a;      ????????????
281
282		UseSweepMethod(d, A, B, C, N, pTable[k]); 
283	}
284}