/work_files/curs/curs/parabolicVector_solver.cpp

http://cm-kp-406.googlecode.com/ · C++ · 198 lines · 157 code · 35 blank · 6 comment · 16 complexity · b231d5e1a0b934191c150004a35f1bcd MD5 · raw file

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