/project/curs/curs/backTime_solver.cpp
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}