/project/curs/curs/reverseTask_solver.cpp

  1#include "reverseTask_solver.h"
  2
  3#include <iomanip>
  4#include <fstream>
  5
  6void swap(double* link1, double* link2)
  7{
  8	double *temp;
  9	temp = link1;
 10	link1 = link2;
 11	link2 = temp;
 12}
 13
 14void vcopy(double* _output, double* _input, int size)
 15{
 16	for (int i = 0; i < size; ++i)
 17	{
 18		_output[i] = _input[i];
 19	}
 20}
 21
 22
 23CReverseTask_solver::CReverseTask_solver()
 24{
 25	iter_num = 0;
 26	cur_q = NULL;
 27	prev_q = NULL;
 28	curGr = NULL;
 29	prevGr = NULL;
 30	curS = NULL;
 31	prevS = NULL;
 32	accuracy = 0.01; 
 33}
 34
 35CReverseTask_solver::~CReverseTask_solver()
 36{
 37	if (cur_q != NULL)
 38		delete [] cur_q;
 39}
 40
 41void CReverseTask_solver::solve_task()
 42{
 43	stepT = maxT / K;
 44
 45	//??? phi1 ??????????? ???????? q ? ?? ?????? ?????!! (q - ??????? ???????? ???????)
 46	CParabolicSolver* direct_solver = new CParabolicSolver();
 47	direct_solver->setGrid(N, K, xb, xe, maxT);
 48	direct_solver->setTask(a, b, c, f, alpha1, betta1, alpha2, betta2, phi1, phi2, psi);
 49	direct_solver->solve_implicit();
 50	direct_solver->CheckSolve("direct_solve.txt");
 51	directTable = direct_solver->GetTable();
 52	analytic = direct_solver->getLeftBorder();
 53
 54	CParabolicVector_solver* directVector_solver = new CParabolicVector_solver();
 55	directVAq_solver = new CParabolicVector_solver();
 56	directVAs_solver = new CParabolicVector_solver();
 57	CBackTimeSolve* backTime_solver = new CBackTimeSolve();
 58
 59	cur_q = new double[K+1];
 60	prev_q = new double[K+1];
 61	curGr = new double[K+1];
 62	prevGr = new double[K+1];
 63	curS = new double[K+1];
 64	prevS = new double[K+1];
 65
 66	// ????????? ??????????? ??????? q
 67	for (int i = 0; i < K + 1; ++i)
 68	{
 69		cur_q[i] = 1.0;
 70	}
 71
 72	// ?????? ??? ?????
 73		//????????? cur_q 
 74
 75	// ?????? ??? ???
 76
 77	//?????????? ?????????
 78
 79	//curGr = ... ???????? ????????
 80
 81	directVector_solver->setGrid(N, K, xb, xe, maxT);
 82	//????????? ????????? (????? ?????...)
 83	directVector_solver->setTask(a, b, c, f, task_lambda, 0, 0, 1, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????
 84//	directVector_solver->setTask(a, b, c, f, task_lambda, 0, 1, 0, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????
 85	directVector_solver->solve_implicit();
 86
 87	print_table(directVector_solver->getTable(), "directVectorGradient.txt");
 88
 89	double* right_border = directVector_solver->getRightBorder();
 90
 91	//????????? ????????? (????? ?????...) ???? ??????????? ??? ??????? ? ?????? ???????.
 92	backTime_solver->setTask(a, b, c, task_lambda, 0, task_lambda, 0, xb, xe, maxT, N, K, right_border, analytic);
 93	backTime_solver->getSolve();
 94	//print_table(backTime_solver->getTable(), "backTimeGradient.txt");
 95	backTime_solver->getGradient(curGr);
 96
 97	iter_gamma = 0.0;
 98
 99	vcopy(curS, curGr, K + 1);	//curS = curGr;	// ????????? ???????????? //??????????? ???????????! ??? ???????????????? ??????!
100
101	iter_betta = get_betta(curGr, curS);
102
103	vcopy(prev_q, cur_q, K + 1);		// ???????? ??????????
104	get_q(curS, prev_q, iter_betta, cur_q);
105	//????? ??? k = 0: curS_0, cur_q_1, prev_q_0, curGr_0
106
107	double eps;
108	double prev_eps = 1.0;
109	eps = get_discrepancy(cur_q, analytic);
110
111	while ( eps > accuracy )
112	{
113		iter_num++;
114	
115		vcopy(prevS, curS, K + 1);		//prevS = curS;			//???????? ?????????? S		// ??????
116		vcopy(prevGr, curGr, K + 1);	//prevGr = curGr;		//???????? ?????????? ????????
117
118		// ?????? ??? ?????
119		//????????? cur_q 
120
121		// ?????? ??? ???
122
123		//?????????? ?????????
124
125		//curGr = ... ???????? ????????
126		directVector_solver->setGrid(N, K, xb, xe, maxT);
127		//????????? ????????? (????? ?????...)
128//		directVector_solver->setTask(a, b, c, f, task_lambda, 0, 1, 0, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????
129		directVector_solver->setTask(a, b, c, f, task_lambda, 0, 0, 1, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????
130		directVector_solver->solve_implicit();
131		double* right_border = directVector_solver->getRightBorder();
132
133		//????????? ????????? (????? ?????...) ???? ??????????? ??? ??????? ? ?????? ???????.
134		backTime_solver->setTask(a, b, c, task_lambda, 0, task_lambda, 0, xb, xe, maxT, N, K, right_border, analytic);
135		backTime_solver->getSolve();
136		backTime_solver->getGradient(curGr);
137		delete [] right_border;
138
139		iter_gamma = get_gamma(prevGr, curGr);
140		get_S(curGr, prevS, iter_gamma, curS);	//curS - output
141		iter_betta = get_betta(curGr, curS);
142		vcopy(prev_q, cur_q, K + 1);
143		//swap (prev_q, cur_q);		// ???????? ??????????
144		get_q(curS, prev_q, iter_betta, cur_q);		//cur_q - output
145
146		prev_eps = eps;
147		eps = get_discrepancy(cur_q, analytic);
148
149		if (prev_eps < eps)
150		{
151			printf("Bad luck, your process is foolish\n");
152			break;
153		}
154	}
155
156	printf("eps = %lf\n", eps);
157
158	printf("q:\n");
159	for (int k = 0; k < K + 1; ++k)
160		printf("q[%d] = %lf\n", k, cur_q[k]);
161
162	delete [] prev_q;
163	delete [] curGr;
164	delete [] prevGr;
165	delete [] curS;
166	delete [] prevS;
167}
168
169double CReverseTask_solver::get_discrepancy(double* _cur_q, double* _analytic)
170{
171	directVAq_solver->setGrid(N, K, xb, xe, maxT);
172	directVAq_solver->setTask(a, b, c, f, alpha1, betta1, alpha2, betta2, _cur_q, phi2, psi);
173	directVAq_solver->solve_implicit();
174	//print_table(directVAq_solver->getTable(), "q_k+1.txt");
175
176	double* Aq;
177	Aq = directVAq_solver->getLeftBorder();
178	double sum = 0.0;
179
180	for (int i = 0; i < K + 1; ++i)
181	{
182		sum += pow(Aq[i] - _analytic[i], 2.0);
183	}
184
185	sum *= stepT;
186
187	return sum;
188}
189
190double CReverseTask_solver::get_betta(double* _curGr, double* _curS)
191{
192	double answer = 0.0;
193	double* VectorAs;
194
195	directVAs_solver->setGrid(N, K, xb, xe, maxT);
196	directVAs_solver->setTask(a, b, c, f, alpha1, betta1, alpha2, betta2, _curS, phi2, psi);
197	directVAs_solver->solve_implicit();
198	VectorAs = directVAs_solver->getLeftBorder();
199
200	double sum_up = 0.0;
201	double sum_down = 0.0;
202
203	for(int i = 0; i < K + 1; i++)
204	{
205		sum_up += _curGr[i] * _curS[i];
206		sum_down += VectorAs[i] * VectorAs[i];
207	}
208
209	answer = sum_up / (2.0 * sum_down);
210
211	return answer;
212}
213
214double CReverseTask_solver::get_gamma(double* _prevGr, double* _curGr)
215{
216	double ans = 0.0;
217	double sum_cur = 0.0;
218	double sum_prev = 0.0;
219
220	for(int i = 0; i < K + 1; i++)
221	{
222		sum_cur += _curGr[i] * _curGr[i];
223		sum_prev += _prevGr[i] * _prevGr[i];
224	}
225
226	ans = sum_cur / sum_prev;
227	return ans;
228}
229
230void CReverseTask_solver::get_S(double* _curGr, double* _prevS, double _iter_gamma, double* _curS)	// ???????? ?????? ? curS
231{
232	for (int i = 0; i < K + 1; ++i)
233	{
234		_curS[i] = _curGr[i] + _iter_gamma * _prevS[i];
235	}
236}
237
238void CReverseTask_solver::get_q(double* _curS, double* _prev_q, double _iter_betta, double* _cur_q)	// ???????? ?????? ? cur_q
239{
240	for (int i = 0; i < K + 1; ++i)
241	{
242		_cur_q[i] = _prev_q[i] - _iter_betta * _curS[i];
243	}
244}
245
246void CReverseTask_solver::print_table(double** a, string fname)
247{
248	ofstream out(fname.c_str());
249	for (int i = 0; i < K+1; ++i)
250	{
251		for (int j = 0; j < N + 1; ++j)
252		{
253			out<< fixed;
254		out <<  setprecision (5)  <<"t = " << i * stepT <<"   x = " << j * (xe - xb) / N << "   u(x,t) = " << a[i][j] << endl;
255		}
256			out << endl;
257	}
258}
259
260//////////////////////////////////////////////////////////////////////////