API tools faq
paste
Guest User

Untitled

a guest
Jun 5th, 2018
69
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 27.96 KB | None | 0 0
raw download clone embed print report
  1. #include "stdio.h"
  2. #include <stdlib.h>
  3. #include "math.h"
  4. #include "main.h"
  5. // #include "user.c"
  6. // #include "no_lin.c"
  7. // #include "print_f.c"
  8.  
  9. double V_nplus1_m(double v_n_mplus1, double v_n_m);
  10. double V_nplus1_m_nelin(double v_n_mplus1, double v_n_m);
  11. double proverka(double value);
  12. int stepx_by_x(double x);
  13. int stept_by_t(double t);
  14. double t_by_stept(int step);
  15. double x_by_stepx(int step);
  16. void print_block(double * block, int length);
  17. double delta(double v_n_m, double x, double t);
  18. double delta_nel(double v_n_m, double x, double t);
  19. void norma_c(double delta, double * max);
  20. void norma_l(double v_n_m, double *norma);
  21. void print_matrix(double * matrix, int length);
  22. void print_block_to_file(double * block, int length, FILE * fout);
  23. void begin_value(double *block, int length);
  24. void characteristic(double * block);
  25. void characteristic_for_b(double * b);
  26. void next_block(double * block, double step_by_t,double * norma_DELTA_c_line, double * norma_delta_c_line,double * norma_DELTA_l_line, double * norma_delta_l_line);
  27. void next_block_nelin(double * block, double step_by_t,double * norma_DELTA_c_line, double * norma_delta_c_line, double * norma_DELTA_l_line, double * norma_delta_l_line);
  28.  
  29.  
  30. void matrix_b(double * b);
  31. void matrix_b1(double * b);
  32. void matrix_U( double * U_1, double * U_2, double * b);
  33. void matrix_V_hat(double * U_1, double * U_2, double * b, double * V_hat, double step_by_t,//неявный случай
  34. double * norma_DELTA_c_line, double * norma_delta_c_line, //[c]
  35. double * norma_DELTA_l_line, double * norma_delta_l_line);
  36. void equall(double * b, double * V_hat);
  37.  
  38. void equal1(double * b, double * V_hat);
  39. void main_loop_lin_simple(void);
  40. void main_loop_lin_nelin(void);
  41. void begin_value_pristrel(double * block, int length);
  42. void matrix_J( double *U_1, double * U_2, double * V_pr, double *G_2);
  43. void matrix_G(double * V_pr, double * G, double * V);
  44. void matrix_G1(double * V_pr, double * G, double * V);
  45. void new_vector_X(double * U_1, double * U_2, double * G_2, double * X);
  46. void new_V_pr(double * X, double * V_pr);
  47. void equal_G2_G1(double * G_1, double G_2[m-2]);
  48. void equal_V_Vpr(double * V, double * V_pr);
  49. void equal_V_Vpr_print(double * V, double * V_pr, int t);
  50. double norma_G(double * G);
  51. double norma_G_2(double * G);
  52.  
  53. void print_block__xvaluex(double * block, int length, FILE * fout);
  54.  
  55. void equal(double * b, double * V_hat);
  56.  
  57.  
  58. // ----------------------------------------------------------------
  59.  
  60.  
  61. double proverka(double value){
  62. if (fabs(value) < 1e-100) value = 1e-100;
  63. return value;
  64. }
  65.  
  66. int stepx_by_x(double x){
  67. return (int)(fabs(x - x_begin)/h);
  68. }
  69.  
  70. int stept_by_t(double t){
  71. return (int)(fabs(t - t_begin)/tao);
  72. }
  73.  
  74. double t_by_stept(int step){
  75. return (double)(t_begin + step*tao);
  76. }
  77.  
  78. double x_by_stepx(int step){
  79. return (double)(x_begin + step*h);
  80. }
  81.  
  82.  
  83. // ----------------------------------------------------------------
  84.  
  85.  
  86. double V_nplus1_m(double v_n_mplus1, double v_n_m){
  87. return (tao/(2*h))*(v_n_mplus1 - v_n_m) + v_n_m;
  88. }
  89.  
  90. double V_nplus1_m_nelin(double v_n_mplus1, double v_n_m){
  91. return ((tao/(2.0*h))*(v_n_mplus1*v_n_mplus1 - v_n_m*v_n_m) + v_n_m);
  92. }
  93.  
  94. // ----------------------------------------------------------------
  95.  
  96.  
  97. void print_matrix(double * matrix, int length){
  98. int i,j;
  99. printf("\nMATRIX PRINT\n");
  100. for (i = 0; i < length; i++){
  101. for (j = 0; j < length; j++){
  102. printf("%-4.1f ", *(matrix + i*length + j));
  103. }
  104. printf("\n");
  105. }
  106. printf("\n");
  107. }
  108.  
  109. void print_block(double * block, int length){
  110. int i;
  111. for (i = 0; i < length; ++i)
  112. {
  113. printf("%f %-4.5g \n", x_by_stepx(i), block[i]);
  114. }
  115. printf("\n");
  116. }
  117.  
  118. void print_block_to_file(double * block, int length, FILE * fout){
  119. int i;
  120. for (i = 0; i < length; ++i)
  121. {
  122. fprintf(fout, "%-8g ", block[i]);
  123. }
  124. fprintf(fout, "\n");
  125. }
  126.  
  127. void print_block__xvaluex(double * block, int length, FILE * fout){
  128. int i;
  129. for (i = 0; i < length; ++i)
  130. {
  131. fprintf(fout, "%-8g %-8g \n", x_by_stepx(i), block[i]);
  132. }
  133. fprintf(fout, "\n");
  134. }
  135.  
  136.  
  137. // ----------------------------------------------------------------
  138.  
  139.  
  140. double delta(double v_n_m, double x, double t){
  141. double delta;
  142. double eps=1e-6;
  143. if (x<=-0.5*t){
  144. delta = fabs(v_n_m - 1);
  145. // printf("delta %g\n", delta);
  146. }
  147. else if ((fabs(x) < 0.5*t + eps)&&(fabs(x) > 0.5*t - eps)){
  148. delta = fabs(v_n_m - 0.5);
  149. // printf("delta %g\n", delta);
  150. }
  151. else {
  152. delta = fabs(v_n_m - 0);
  153. }
  154. return delta;
  155. }
  156.  
  157. double delta_nel(double v_n_m, double x, double t){
  158. double delta;
  159. // printf("[t] %g\n", t);
  160. if (x<=-t){
  161. delta = fabs(v_n_m - 1);
  162. }
  163. else if((x>-t)&&(x<=0)){
  164. delta = fabs(v_n_m + x/t);
  165. }
  166. else if (x>0) {delta = fabs(v_n_m);
  167. };
  168. return delta;
  169. }
  170.  
  171.  
  172. void norma_c(double delta, double * max){
  173. if (fabs(delta) > *max) *max = fabs(delta);
  174. }
  175.  
  176.  
  177. void norma_l(double v_n_m, double * norma){
  178. *norma += fabs(v_n_m);
  179. }
  180.  
  181. // ----------------------------------------------------------------
  182.  
  183. //double delta_izm (double * izm, double * V_pr, M){
  184. // for
  185. //}
  186.  
  187. // ----------------------------------------------------------------
  188.  
  189. double res(double x, double t){
  190. double res;
  191. double eps = 1e-6;
  192. if (x<=-0.5*t){
  193. res = 1;
  194. }
  195. else if ((fabs(x) < 0.5*t + eps)&&(fabs(x) > 0.5*t - eps)){
  196. res = 0.5;
  197. }
  198. else{
  199. res = 0;
  200. }
  201. return res;
  202. }
  203.  
  204. double res_nel( double x, double t){
  205. double res;
  206. if (x<=-t){
  207. res = 1;
  208. }
  209. else if((x<=0)&&(x>-t)){
  210. res = -x/t;
  211. }
  212. else {res = 0;};
  213.  
  214. return res;
  215. }
  216.  
  217. double res_test(double x, double t){
  218. return (1-x)/(t+2);
  219. }
  220. // ----------------------------------------------------------------
  221.  
  222.  
  223.  
  224. void begin_value(double * block, int length){
  225. int i;
  226. for (i = 0; i < length; ++i)
  227. {
  228. if (x_by_stepx(i) <= 0)
  229. block[i] = 1.0;
  230. else
  231. block[i] = 0.0;
  232. }
  233. }
  234.  
  235.  
  236. void characteristic(double * block){
  237. block[0] = 1.0;
  238. block[m-1] = 0;
  239. }
  240.  
  241.  
  242. void characteristic_for_b(double * b){
  243. b[0] = 1.0 - tao/(4.0*h);
  244. b[m-3] = 0.0;
  245. }
  246.  
  247.  
  248. void next_block(double * block, double step_by_t,
  249. /*явный случай*/
  250. double * norma_DELTA_c_line, double * norma_delta_c_line, //[c]
  251. double * norma_DELTA_l_line, double * norma_delta_l_line //[l]
  252. ){
  253. int i;
  254. double x = x_begin;
  255. double del;
  256. *norma_DELTA_c_line = 0; *norma_delta_c_line = 0; //начальное максимальное значение [c]
  257. *norma_DELTA_l_line = 0; *norma_delta_l_line = 0; //начальное значение [l]
  258.  
  259. characteristic(block);
  260. for (i = 0; i < m-1; ++i)
  261. {
  262. block[i] = V_nplus1_m(block[i+1], block[i]);
  263.  
  264. del = delta(block[i], x, (step_by_t+1)*tao); // разность между численным решением и точным решением
  265. norma_c(del, norma_DELTA_c_line); //[c]
  266. norma_c(block[i], norma_delta_c_line); //[c]
  267. norma_l(del, norma_DELTA_l_line); //[l]
  268. norma_l(block[i], norma_delta_l_line); //[l]
  269.  
  270.  
  271. printf("V: %g x: %g del: %g \n", block[i], x, del);
  272. x+=h;
  273. if (del >= 1) { printf("t: %g\n ", tao*(step_by_t+1));};
  274. }
  275. // printf("NORMA 1 %g\n", *norma_DELTA_c_line);
  276. // printf("NORMA 2 %g\n", *norma_delta_c_line);
  277. }
  278.  
  279.  
  280. void next_block_nelin(double * block, double step_by_t,
  281. /*нелинейный случай*/
  282. double * norma_DELTA_c_line, double * norma_delta_c_line, //[c]
  283. double * norma_DELTA_l_line, double * norma_delta_l_line //[l]
  284. ){
  285. int i;
  286. double x = x_begin;
  287. double del;
  288. *norma_DELTA_c_line = 0; *norma_delta_c_line = 0; //начальное максимальное значение [c]
  289. *norma_DELTA_l_line = 0; *norma_delta_l_line = 0; //начальное значение [l]
  290.  
  291. characteristic(block);
  292. for (i = 0; i < m-1; ++i)
  293. {
  294. block[i] = V_nplus1_m_nelin(block[i+1], block[i]);
  295. del = delta_nel(block[i], x, (step_by_t+1)*tao); // разность между численным решением и точным решением
  296. norma_c(del, norma_DELTA_c_line); //[c]
  297. norma_c(block[i], norma_delta_c_line); //[c]
  298. norma_l(del, norma_DELTA_l_line); //[l]
  299. norma_l(block[i], norma_delta_l_line); //[l]
  300. // printf("V: %g x: %g del: %g \n", block[i], x, del);
  301. printf("%g %g \n", x, del);
  302. if (del >= 1) { printf("t: %g\n ", tao*(step_by_t+1));};
  303. x+=h;
  304. }
  305. printf("kek\n");
  306. // printf("NORMA 1 %g\n", *norma_DELTA_c_line);
  307. // printf("NORMA 2 %g\n", *norma_delta_c_line);
  308.  
  309. }
  310.  
  311. // ----------------------------------------------------------------
  312.  
  313.  
  314. void matrix_b(double * b){
  315. int i;
  316. characteristic_for_b(b);
  317. for (i = 1; i < m-1; i++)
  318. {
  319. if (x_by_stepx(i) <= 0)
  320. b[i] = 1.0;
  321. else
  322. b[i] = 0.0;
  323. }
  324. }
  325.  
  326. void matrix_b1(double * b){
  327. /* omega = 1 */
  328. int i;
  329. characteristic_for_b(b);
  330. for (i = 1; i < m-2; i++)
  331. {
  332. if ((x_by_stepx(i+2)==0)&&(x_by_stepx(i+1)!=0)){ // ÎÁßÇÀÒÅËÜÍÎÅ ÏÎÏÀÄÀÍÈÅ Â 0
  333. b[i] = 1 - tao/8;
  334. } else if ((x_by_stepx(i+1)==0)&&(x_by_stepx(i)!=0)){
  335. b[i] = 1 + (3*tao)/8;
  336. } else if ((x_by_stepx(i)==0)&&(x_by_stepx(i-1)!=0)){
  337. b[i] = -(3*tao)/8;
  338. } else if ((x_by_stepx(i-1)==0)&&(x_by_stepx(i-2)!=0)){
  339. b[i] = tao/8;
  340. } else if (x_by_stepx(i+2) < 0){
  341. b[i] = 1.0;
  342. } else
  343. b[i] = 0.0;
  344. }
  345. }
  346.  
  347. // ----------------------------------------------------------------
  348.  
  349. void matrix_U( double * U_1, double * U_2, double * b){
  350. int i;
  351. double l=0; // íà÷àëüíîå l
  352. for (i=0; i<m-2; i++){
  353. if (i==0) { //âûðîæäåííûé ñëó÷àé
  354. U_1[i] = 1;
  355. U_2[i] = -tao/(4*h);
  356. b[i] = b[i];
  357. l = (tao/(4*h))/U_1[i];
  358. } else if (i!=m-3) {
  359. U_1[i] = 1 - l*U_2[i-1];
  360. U_2[i] = -tao/(4*h);
  361. b[i] = b[i] - l*b[i-1];
  362. l = (tao/(4*h))/U_1[i];
  363. } else {
  364. U_1[i] = 1 - l*U_2[i-1];
  365. b[i] = b[i] - l*b[i-1];
  366. }
  367.  
  368. }
  369. }
  370.  
  371. void matrix_V_hat(double * U_1, double * U_2, double * b, double * V_hat, double step_by_t, //неявный случай
  372. double * norma_DELTA_c_line, double * norma_delta_c_line,
  373. double * norma_DELTA_l_line, double * norma_delta_l_line){
  374.  
  375. int i;
  376. double del;
  377. *norma_DELTA_c_line = 0; *norma_delta_c_line = 0; //начальное максимальное значение [c]
  378. *norma_DELTA_l_line = 0; *norma_delta_l_line = 0;
  379.  
  380.  
  381. characteristic(V_hat);
  382. //printf("V_hat[%d]: %g\n", m-1,V_hat[m-1]);
  383. for(i=m-2; i>=1; i--){
  384. if (i==m-2) V_hat[i] = b[i-1]/U_1[i-1];
  385. else {
  386. V_hat[i] = (b[i-1] - U_2[i-1]*V_hat[i+1])/U_1[i-1];
  387. }
  388.  
  389. del = delta(V_hat[i], x_by_stepx(i), (step_by_t+1)*tao); // разность между численным решением и точным решением
  390. norma_c(del, norma_DELTA_c_line); //[c]
  391. norma_c(V_hat[i], norma_delta_c_line); //[c]
  392. norma_l(del, norma_DELTA_l_line);
  393. norma_l(V_hat[i], norma_delta_l_line);
  394. printf("%g %g \n", x_by_stepx(i), del);
  395. // printf("NORMA_DELTA_c_line %g\n", *norma_DELTA_c_line);
  396. // printf("V_hat[%d]: %g\n", i,V_hat[i]);
  397. }
  398. }
  399.  
  400. // ----------------------------------------------------------------
  401.  
  402.  
  403. void equal(double * b, double * V_hat){
  404. /* omega = 0 */
  405. int i;
  406. b[0] = V_hat[1] - tao/(4*h);
  407. for(i=1; i<m-2; i++){
  408. b[i] = V_hat[i+1];
  409. }
  410. }
  411.  
  412. void equal1(double * b, double * V_hat){
  413. /* omega = 1 */
  414. int i;
  415. b[0] = V_hat[1] - (tao/8.)*(V_hat[3] - 4*V_hat[2] + 6*V_hat[1] - 3.0) - tao/(4*h);
  416. b[1] = V_hat[2] - (tao/8.)*(V_hat[4] - 4*V_hat[3] + 6*V_hat[2] - 4.0*V_hat[1] + 1);
  417. for(i=2; i<m-4; i++){
  418. b[i] = V_hat[i+1] - (tao/8.)*(V_hat[i+3] - 4*V_hat[i+2] + 6*V_hat[i+1] - 4*V_hat[i] + V_hat[i-1]);
  419. }
  420. b[m-4] = V_hat[m-3] - (tao/8.)*(V_hat[m-5] - 4*V_hat[m-4] + 6*V_hat[m-3] - 4*V_hat[m-2]);
  421. b[m-3] = V_hat[m-2] - (tao/8.)*(V_hat[m-4] - 4*V_hat[m-3] + 6*V_hat[m-2]);
  422. }
  423.  
  424. /* ---------------- Явные схемы --------------------------------- */
  425.  
  426. void main_loop_lin_simple(){
  427.  
  428. int i;
  429. double block[m];
  430. double norma_DELTA_c_line; double norma_delta_c_line;
  431. double norma_DELTA_l_line; double norma_delta_l_line;
  432. FILE * fout = fopen("output_1.txt", "w+");
  433.  
  434. int length = (int)(sizeof(block)/sizeof(double));
  435. begin_value(block, length);
  436.  
  437. for (i = 0; i < n; ++i){
  438. next_block(block, i, &norma_DELTA_c_line, &norma_delta_c_line, &norma_DELTA_l_line, &norma_delta_l_line);
  439. /*norma_c(norma_DELTA_c_line, &norma_DELTA_c);
  440. norma_c(norma_delta_c_line, &norma_delta_c);
  441. norma_l(norma_DELTA_l_line, &norma_DELTA_l);
  442. norma_l(norma_delta_l_line, &norma_delta_l);*/
  443.  
  444. print_block_to_file(block, length, fout);
  445.  
  446. }
  447.  
  448. printf("NORMA_DELTA_c %g\n", norma_DELTA_c_line);
  449. printf("NORMA_delta_c %g\n", norma_DELTA_c_line/norma_delta_c_line);
  450. printf("NORMA_DELTA_l %g\n", h*norma_DELTA_l_line);
  451. printf("NORMA_delta_l %g\n", norma_DELTA_l_line/norma_delta_l_line);
  452.  
  453. printf("LIN__________\n");
  454. printf("tao h DeltaC DeltaL deltaC delta l\n");
  455. printf("%g & %g & %e & %e & %e & %e \n", tao, h, norma_DELTA_c_line, h*norma_DELTA_l_line, norma_DELTA_c_line/norma_delta_c_line, norma_DELTA_l_line/norma_delta_l_line);
  456.  
  457. }
  458.  
  459. /* ---------------- Нелинейный случай (неявный) -------------------*/
  460.  
  461. void main_loop_lin_nelin(){
  462. int i;
  463. double block[m];
  464. double norma_DELTA_c_line; double norma_delta_c_line;
  465. double norma_DELTA_l_line; double norma_delta_l_line;
  466. FILE * fout = fopen("output_1.txt", "w+");
  467.  
  468. int length = (int)(sizeof(block)/sizeof(double));
  469. begin_value(block, length);
  470.  
  471. for (i = 0; i < n; ++i){
  472. next_block_nelin(block, i, &norma_DELTA_c_line, &norma_delta_c_line, &norma_DELTA_l_line, &norma_delta_l_line);
  473.  
  474. // print_block_to_file(block, length, fout);
  475.  
  476. }
  477.  
  478. print_block__xvaluex(block, length, fout);
  479.  
  480.  
  481.  
  482. printf("NORMA_DELTA_c %g\n", norma_DELTA_c_line);
  483. printf("NORMA_delta_c %g\n", norma_DELTA_c_line/norma_delta_c_line);
  484. printf("NORMA_DELTA_l %g\n", h*norma_DELTA_l_line);
  485. printf("NORMA_delta_l %g\n", norma_DELTA_l_line/norma_delta_l_line);
  486.  
  487.  
  488. printf("NO LIN__________\n");
  489. printf("tao%-7s h%-9s %-7sDeltaC %-7sDeltaL %-7sdeltaC %-7sdelta l\n", "","","","","","");
  490. printf("%-10g & %-10g & %-10e & %-10e & %-10e & %-10e \n\n", tao, h, norma_DELTA_c_line, h*norma_DELTA_l_line, norma_DELTA_c_line/norma_delta_c_line, norma_DELTA_l_line/norma_delta_l_line);
  491.  
  492. }
  493.  
  494. /* ---------------- Нелинейный случай (неявный) -------------------*/
  495.  
  496. void begin_value_pristrel(double * block, int length){
  497. int i;
  498. for (i = 0; i < length; ++i){
  499. block[i] = 1.0; //-x_by_stepx(i)*0.5+0.5;
  500. }
  501. }
  502.  
  503. void matrix_J( double * U_1, double * U_2, double * V_pr, double * G_2){
  504. int i;
  505. double l = 0; // íà÷àëüíîå l
  506. double a = 2*tao/(4*h);
  507.  
  508. for (i=0; i<m-2; i++){
  509. if (i==0) { //âûðîæäåííûé ñëó÷àé
  510. U_1[i] = 1;
  511. U_2[i] = proverka(-a*V_pr[2]);
  512. G_2[i] = proverka(G_2[i]);
  513. l = proverka(a*V_pr[1]/U_1[i]); // ???????? 1==0?
  514. } else if (i!=m-3) {
  515. U_1[i] = proverka(1 - l*U_2[i-1]);
  516. U_2[i] = proverka(-a*V_pr[i+2]);
  517. G_2[i] = proverka(G_2[i] - l*G_2[i-1]);
  518. l = proverka(a*V_pr[i+1]/U_1[i]); // ??????????
  519. } else {
  520. U_1[i] = proverka(1 - l*U_2[i-1]);
  521. G_2[i] = proverka(G_2[i] - l*G_2[i-1]);
  522. }
  523. }
  524.  
  525.  
  526.  
  527. // for (int i = 0; i < m-3; ++i){
  528. // printf("U_1[%-2d]: %-13g, U_2[%-2d]: %-13g \n", i, U_1[i], i, U_2[i]);
  529. // }
  530. }
  531.  
  532. void matrix_G(double * V_pr, double * G, double * V){
  533. int i;
  534. double a = tao/(4*h);
  535.  
  536. for (i=0; i<m-2; i++){
  537. if (i==0) {
  538. G[0] = -(V_pr[1] - V[1] - a*(V_pr[2]*V_pr[2] - 1));
  539. } else if (i!=m-3) {
  540. G[i] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i]));
  541. } else {
  542. G[m-3] = -(V_pr[m-2] - V[m-2] - a*(0 - V_pr[m-3]*V_pr[m-3]));
  543. }
  544.  
  545. }
  546. }
  547.  
  548. void matrix_G_w1(double * V_pr, double * G, double * V){
  549. int i;
  550. double a = tao/(4*h); //111111111111111111111100000000000000000000000
  551.  
  552. for (i=0; i<m-2; i++){
  553. if (i==0) {
  554. G[0] = -(V_pr[1] - V[1] - a*(V_pr[2]*V_pr[2] - 1) - tao*0.125*(-3 + V[3] - 4*V[2] + 6*V[1] )) ;
  555. } else if (i==1){
  556. G[1] = -(V_pr[2] - V[2] - a*(V_pr[3]*V_pr[3] - V_pr[1]*V_pr[1]) - tao*0.125*(V[i+3] - 4*V[i+2] + 6*V[i+1] - 4*V[i] + 1));
  557. }
  558. else if ((i!=m-3)&&(i!=m-4)) {
  559. G[i] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i]) - tao*0.125*(V[i+3] - 4*V[i+2] + 6*V[i+1] - 4*V[i] + V[i-1]));
  560. } else if (i==m-4){
  561. G[m-4] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i]) - tao*0.125*(0 - 4*V[i+2] + 6*V[i+1] - 4*V[i] + V[i-1]));
  562. } else if (i==m-3){
  563. G[m-3] = -(V_pr[i+1] - V[i+1] - a*(0 - V_pr[i]*V_pr[i]) - tao*0.125*(0 - 0 + 6*V[i+1] - 4*V[i] + V[i-1]));
  564. }
  565.  
  566. }
  567. }
  568.  
  569. void new_vector_X(double * U_1, double * U_2, double * b, double * X){
  570.  
  571. int i;
  572.  
  573. for(i=m-3; i>=0; i--){
  574. if (i == m-3) {
  575.  
  576. X[i] = proverka(b[i]/U_1[i]);
  577. }
  578. else {
  579. X[i] = proverka((b[i] - U_2[i]*X[i+1])/U_1[i]);
  580. }
  581. }
  582. }
  583.  
  584. void new_V_pr(double * X, double * V_pr){
  585. for (int i = 0; i < m-2; ++i){
  586. V_pr[i+1] = V_pr[i+1] + X[i];
  587. }
  588. }
  589.  
  590. void equal_G2_G1(double * G_1, double * G_2){
  591. for (int i = 0; i < m-2; ++i)
  592. {
  593. G_2[i] = G_1[i];
  594. }
  595. }
  596.  
  597. void equal_V_Vpr(double * V, double * V_pr){
  598. for (int i = 0; i < m; ++i)
  599. {
  600. V[i] = V_pr[i];
  601. }
  602. }
  603.  
  604. double norma_G(double * G){
  605. double norma = 0;
  606. for (int i = 0; i < m-2; ++i){
  607. norma += fabs(G[i]);
  608. }
  609.  
  610. return norma*h;
  611. }
  612.  
  613. /* ------------------ ТЕСТ -----------------*/
  614.  
  615. double delta_test(double v_n_m, double x, double t){
  616. double delta;
  617. delta = fabs(v_n_m-((1- x)/(t+2.0)));
  618. return delta;
  619. }
  620.  
  621. void begin_value_test(double block[m], int length){
  622. int i;
  623. for (i = 0; i < length; ++i)
  624. {
  625. block[i] = -0.5*x_by_stepx(i)+0.5;
  626. }
  627. //printf("BLOCK[%g] = %g\n", x_by_stepx(i), block[i]);
  628.  
  629.  
  630. }
  631.  
  632. void characteristic_test(double block[m], double t){
  633. block[0] = 2.0/(2.0+t);
  634. block[m-1] = 0.0;
  635. }
  636.  
  637. void matrix_G_test(double V_pr[m], double G[m-2], double V[m], double t){
  638. int i;
  639. double a = tao/(4.0*h);
  640.  
  641. for (i=0; i<m-2; i++){
  642. if (i==0) {
  643. G[0] = -(V_pr[1] - V[1] - a*(V_pr[2]*V_pr[2] - (2.0/(2.0+t))*(2.0/(2.0+t))));
  644. } else if (i==m-3) {
  645. G[m-3] = -(V_pr[m-2] - V[m-2] - a*(0.0- V_pr[m-3]*V_pr[m-3]));
  646. } else {G[i] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i]));
  647.  
  648. }
  649.  
  650. }
  651. }
  652.  
  653. void matrix_G1_test(double V_pr[m], double G[m-2], double V[m], double t){
  654. int i;
  655. double a = tao/(4.0*h);
  656.  
  657. for (i=0; i<m-2; i++){
  658. if (i==0) {
  659. G[0] = -(V_pr[1] - V[1] - a*(V_pr[2]*V_pr[2] - (2.0/(2.0+t))*(2.0/(2.0+t))) -tao*0.125*(V[3]-4*V[2]+6*V[1]-4*(2.0/(2.0+t))+(3+h)/(t+2)));
  660. } else if (i == 1)
  661. {
  662. G[1] = -(V_pr[2] - V[2] - a*(V_pr[3]*V_pr[3] - V_pr[1]*V_pr[1]) -tao*0.125*(V[4]-4*V[3]+6*V[2]-4*V[1]+(2.0/(2.0+t))));
  663. }
  664. else if ((i!=m-3)&&(i!=m-4) ){
  665. G[i] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i]) - tao*0.125*(V[i+3]-4*V[i+2]+6*V[i+1]-4*V[i]+V[i-1]));
  666. }
  667. else if (i== m-4){
  668. G[i] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i])-tao*0.125*(0-4*V[i+3]+6*V[i+1]-4*V[i]+V[i-1]));
  669. }
  670. else { //i = m-3
  671. G[m-3] = -(V_pr[m-2] - V[m-2] - a*(0- V_pr[m-3]*V_pr[m-3])- tao*0.125*((-h)/(t+2)+6*V[m-2]-4*V[m-3]+V[m-4]));
  672. }
  673.  
  674. }
  675.  
  676. }
  677.  
  678. void equal_V_Vpr_print(double V[m], double V_pr[m], int t){
  679. int i;
  680. double del;
  681. for (i = 0; i < m-1; ++i)
  682. {
  683. V[i] = V_pr[i];
  684. del = delta_test(V[i], x_by_stepx(i), t);
  685. if(t>0.9) printf("%g %g %g\n", V[i], del, x_by_stepx(i));
  686. }
  687. }
  688.  
  689. void test(void){
  690.  
  691. int i;
  692.  
  693. double V_hat[m];
  694. double U_1[m-2]; double U_2[m-3];
  695.  
  696. int length = (int)(sizeof(V_hat)/sizeof(double));
  697.  
  698. //по всей сетке [c]
  699. double norma_DELTA_c_line; double norma_delta_c_line; //по линиям
  700. //по всей сетке [l]
  701. double norma_DELTA_l_line; double norma_delta_l_line; //по линиям
  702.  
  703. FILE * fout_test = fopen("test.txt", "w");
  704.  
  705. //################ ТЕСТ ##########################
  706.  
  707. int j; double del;double V_pr[m]; double V_0[m]; double G_1[m-2]; double G_2[m-2]; double V[m]; double X[m-2]; double norma_for_G = 1;
  708.  
  709. begin_value_test(V, length);
  710. equal_V_Vpr(V_0, V);
  711.  
  712. for (i = 0; i < n; ++i)
  713. {
  714.  
  715.  
  716. equal_V_Vpr(V_pr, V_0);
  717. characteristic_test(V_pr,(i+1)*tao);
  718.  
  719. while (fabs(norma_for_G) > 0.00001) {
  720.  
  721. printf("i %d\n", i);
  722. printf("%e - NORMAAAAAAAAAAAAAAAAAAAAAAAAAA GGGGGGGGGGGGG\n\n", norma_for_G);
  723.  
  724. matrix_G_test(V_pr, G_1, V, (i+1)*tao );
  725.  
  726.  
  727.  
  728. equal_G2_G1(G_1, G_2);
  729.  
  730. matrix_J(U_1, U_2, V_pr, G_2);
  731.  
  732.  
  733. new_vector_X(U_1, U_2, G_2, X);
  734. new_V_pr(X, V_pr);
  735.  
  736. matrix_G_test(V_pr, G_1, V, (i+1)*tao);
  737. norma_for_G = norma_G(G_1);
  738.  
  739. }
  740.  
  741. norma_for_G = 1;
  742. equal_V_Vpr(V_0,V);
  743. equal_V_Vpr_print(V, V_pr, (i+1)*tao);
  744.  
  745. printf("------------------------------`");
  746.  
  747.  
  748. norma_DELTA_c_line = 0; norma_delta_c_line = 0; //начальное максимальное значение [c]
  749. norma_DELTA_l_line = 0; norma_delta_l_line = 0; //начальное значение [l]
  750.  
  751. for (j = 0; j < m-1; ++j)
  752. {
  753. del = delta_test(V_pr[j], x_by_stepx(j), (i+1)*tao); // разность между численным решением и точным решением
  754.  
  755. norma_c(del, &norma_DELTA_c_line); //[c]
  756. norma_c(V_pr[j], &norma_delta_c_line); //[c]
  757. norma_l(del, &norma_DELTA_l_line); //[l]
  758. norma_l(V_pr[j], &norma_delta_l_line); //[l]
  759. }
  760. }
  761.  
  762. print_block__xvaluex(V_pr, m, fout_test);
  763. printf("tao h DeltaC DeltaL deltaC delta l\n");
  764. printf("%g & %g & %e & %e & %e & %e \n", tao, h, norma_DELTA_c_line, h*norma_DELTA_l_line, norma_DELTA_c_line/norma_delta_c_line, norma_DELTA_l_line/norma_delta_l_line);
  765.  
  766.  
  767.  
  768.  
  769. }
  770.  
  771. int main(){
  772.  
  773.  
  774. // INIT FOR IMPLICIT CASE
  775.  
  776. double U_1[m-2]; double U_2[m-3]; double b[m-2]; double V_hat[m];
  777. double V_pr[m]; double G_1[m-2]; double G_2[m-2]; double V[m]; double X[m-2]; double norma_for_G = 1;
  778. double tochno_lin[m];
  779. double tochno_nelin[m];
  780. double tochno_test[m];
  781. double izm[m];
  782.  
  783. int length = (int)(sizeof(V_hat)/sizeof(double)); // equal m
  784.  
  785. double norma_DELTA_c_line; double norma_delta_c_line;
  786. double norma_DELTA_l_line; double norma_delta_l_line;
  787.  
  788. FILE * fout = fopen("output_2.txt", "w+");
  789. FILE * fout_once = fopen("output_once.txt", "w");
  790. FILE * f_lin = fopen("lin.txt", "w");
  791. FILE * f_nelin = fopen("nelin.txt", "w");
  792. FILE * f_test = fopen("test.txt", "w");
  793. FILE * f_izm = fopen("izm.txt", "w");
  794.  
  795. double del; // нужно для нормы
  796.  
  797. printf("Size: %d x %d\n", m,n);
  798.  
  799. for(int i; i<m; i++){
  800. tochno_lin[i] = res(x_by_stepx(i), 1);
  801. fprintf(f_lin, "%g %g\n", x_by_stepx(i), tochno_lin[i]);
  802. }
  803.  
  804. for(int i; i<m; i++){
  805. tochno_nelin[i] = res_nel(x_by_stepx(i), 1);
  806. fprintf(f_nelin, "%g %g\n", x_by_stepx(i), tochno_nelin[i]);
  807. }
  808.  
  809. for(int i; i<m; i++){
  810. tochno_test[i] = res_test(x_by_stepx(i), 1);
  811. fprintf(f_test, "%g %g\n", x_by_stepx(i), tochno_test[i]);
  812. }
  813.  
  814. /*#########################
  815. Тест
  816. ###########################*/
  817.  
  818. test();
  819.  
  820. //return 0;
  821.  
  822. /*#########################
  823. Явная схема
  824. ###########################*/
  825.  
  826.  
  827. //main_loop_lin_simple();
  828. //main_loop_lin_nelin();
  829. //return 0;
  830.  
  831. /*#########################
  832. Линейный случай
  833. ###########################*/
  834.  
  835. /*
  836. matrix_b1(b);
  837.  
  838. for(int i=0; i<n; i++){
  839. //printf("\nTao: %f, Step %d\n\n", t_by_stept(i+1), i);
  840.  
  841. matrix_U(U_1, U_2, b);
  842. matrix_V_hat(U_1, U_2, b, V_hat, i, &norma_DELTA_c_line, &norma_delta_c_line,&norma_DELTA_l_line, &norma_delta_l_line);
  843. equal1(b, V_hat);
  844. print_block_to_file(V_hat, length, fout);
  845.  
  846. //printf("NORMA_DELTA_c_line %g\n", norma_DELTA_c_line);
  847. //printf("NORMA_DELTA_c %g\n", norma_DELTA_c);
  848.  
  849. }
  850.  
  851. print_block__xvaluex(V_hat, length, fout_once);
  852.  
  853. //printf("NORMA_DELTA_c %g\n", norma_DELTA_c);
  854. //printf("NORMA_delta_c %g\n", norma_DELTA_c/norma_delta_c);
  855. //printf("NORMA_DELTA_l %g\n", h*norma_DELTA_l);
  856.  
  857. printf("NELIN \n");
  858. printf("tao h DeltaC DeltaL deltaC delta l\n");
  859. printf("%g & %g & %e & %e & %e & %e \n", tao, h, norma_DELTA_c_line, h*norma_DELTA_l_line, norma_DELTA_c_line/norma_delta_c_line, norma_DELTA_l_line/norma_delta_l_line);
  860.  
  861. return 0;
  862. */
  863.  
  864. /*#########################
  865. Нелинейный случай
  866. ###########################*/
  867.  
  868. for(int k=0; k<4; k++){
  869.  
  870.  
  871.  
  872. begin_value(V, length);
  873.  
  874. for (int i = 0; i < n; ++i)
  875. {
  876.  
  877. begin_value_pristrel(V_pr, length);
  878. V_pr[0] = 1; //1/(1 + 0.5*(tao*(i+1)))
  879. V_pr[m-1] = 0;
  880.  
  881. while (norma_for_G > 1e-13) {
  882.  
  883. // printf("%g - NORMAAAAAAAAAAAAAAAAAAAAAAAAAA GGGGGGGGGGGGG\n\n", norma_for_G);
  884.  
  885. matrix_G(V_pr, G_1, V);
  886. equal_G2_G1(G_1, G_2);
  887. matrix_J(U_1, U_2, V_pr, G_2);
  888. new_vector_X(U_1, U_2, G_2, X);
  889. new_V_pr(X, V_pr);
  890. matrix_G(V_pr, G_1, V); // изменение G для нового вектора пристреленного
  891.  
  892. // printf("%s", "lol ");
  893. // printf("%s", "kek ");
  894. norma_for_G = norma_G(G_1);
  895. printf("%g\n", norma_for_G);
  896. }
  897.  
  898. norma_for_G = 1;
  899.  
  900. }
  901.  
  902. norma_DELTA_c_line = 0; norma_delta_c_line = 0; //начальное максимальное значение [c]
  903. norma_DELTA_l_line = 0; norma_delta_l_line = 0; //начальное значение [l]
  904.  
  905. for (int j = 0; j < m-1; ++j)
  906. {
  907. del = delta_nel(V_pr[j], x_by_stepx(j), (i+1)*tao); // разность между численным решением и точным решением
  908. printf("%s %f %f %f \n", "X and T and del",x_by_stepx(j), (i+1)*tao, del);
  909. norma_c(del, &norma_DELTA_c_line); //[c]
  910. norma_c(V_pr[j], &norma_delta_c_line); //[c]
  911. norma_l(del, &norma_DELTA_l_line); //[l]
  912. norma_l(V_pr[j], &norma_delta_l_line); //[l]
  913. }
  914.  
  915. print_block(V_pr, m);
  916. equal_V_Vpr(V, V_pr);
  917. printf("I : %d\n", i);
  918.  
  919. equal_V_Vpr(izm, V_pr);
  920. if (k>=1){
  921.  
  922. }
  923. k++;
  924.  
  925. print_block__xvaluex(V_pr, m, fout_once);
  926.  
  927. printf("NELIN \n");
  928. printf("tao h DeltaC DeltaL deltaC delta l\n");
  929. printf("%g & %g & %e & %e & %e & %e \n", tao, h, norma_DELTA_c_line, h*norma_DELTA_l_line, norma_DELTA_c_line/norma_delta_c_line, norma_DELTA_l_line/norma_delta_l_line);
  930. }
  931.  
  932. return 0;
  933. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!