API tools faq
paste
Advertisement
Guest User

Untitled

a guest
May 25th, 2018
65
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 20.75 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.  
  6. double V_nplus1_m(double v_n_mplus1, double v_n_m){
  7. return (tao/(2*h))*(v_n_mplus1 - v_n_m) + v_n_m;
  8. }
  9.  
  10. double V_nplus1_m_nelin(double v_n_mplus1, double v_n_m){
  11. // if (t_local > 0.995)
  12. //{ printf(" \n res2: %.25lf \n v_n_m: %.25lf \n res %.25lf \n \n",(tao*0.5)*(pow(v_n_mplus1, 2) - pow(v_n_m, 2))/h, v_n_m, (tao/(2*h))*(v_n_mplus1*v_n_mplus1 - v_n_m*v_n_m) + v_n_m);};
  13. return ((tao*0.5)*(pow(v_n_mplus1, 2) - pow(v_n_m, 2))/h + v_n_m);
  14. }
  15.  
  16.  
  17. double proverka(double value){
  18. if (fabs(value)<1e-15) value = 0.0;
  19. return value;
  20. }
  21.  
  22. int stepx_by_x(double x){
  23. return (int)(fabs(x - x_begin)/h);
  24. }
  25.  
  26. int stept_by_t(double t){
  27. return (int)(fabs(t - t_begin)/tao);
  28. }
  29.  
  30. double t_by_stept(int step){
  31. return (double)(t_begin + step*tao);
  32. }
  33.  
  34. double x_by_stepx(int step){
  35.  
  36. return (double)((x_begin + step*h));
  37. }
  38.  
  39. void print_block(double * block, int length){
  40. int i;
  41. for (i = 0; i < length; ++i)
  42. {
  43. printf("%-4.5g %g \n", block[i], x_by_stepx(i));
  44. }
  45. printf("\n");
  46. }
  47.  
  48.  
  49. double delta(double v_n_m, double x, double t){
  50. double delta;
  51. /*if (x<=-0.5*t){
  52. delta = fabs(v_n_m - 1);
  53. // printf("delta %g\n", delta);
  54. }
  55. else{
  56. delta = fabs(v_n_m);
  57. // printf("delta %g\n", delta);
  58. }*/
  59.  
  60. /*if (x<=-0.5*t){
  61. delta = fabs(v_n_m);
  62. // printf("delta %g\n", delta);
  63. }
  64. else if ((x>-0.5*t)&&(x<=-0.5*t+0.25)) {
  65. delta = fabs(v_n_m - 4*x-2*t);
  66. // printf("delta %g\n", delta);
  67. }
  68. else {delta = fabs(v_n_m-1);}
  69. return delta;*/
  70.  
  71. if (x<=-0.5*t){
  72. delta = fabs(v_n_m);
  73. // printf("delta %g\n", delta);
  74. }
  75. else if ((x>-0.5*t)&&(x<=-0.5*t+0.25)) {
  76. delta = fabs(v_n_m - 4*x-2*t);
  77. // printf("delta %g\n", delta);
  78. }
  79. else {delta = fabs(v_n_m-1);}
  80. return delta;
  81. }
  82.  
  83. double delta_nel(double v_n_m, double x, double t){
  84. double delta;
  85. /*if (x<=-t){
  86. delta = fabs(v_n_m - 1);
  87. // printf("delta %g\n", delta);
  88. }
  89. else if((x<=0)&&(x>-t)){
  90. delta = fabs(v_n_m+x/t);
  91. // printf("delta %g\n", delta);
  92. }
  93. else {delta = fabs(v_n_m);};*/
  94.  
  95. if ((t<=0.25)&&(x<=0)){
  96. delta = fabs(v_n_m);
  97. // printf("delta %g\n", delta);
  98. }
  99. else if((t<=0.25)&&(x>0)&&(x<=-t+0.25)){
  100. delta = fabs(v_n_m-(4.0*x/(1-4.0*t)));
  101. // printf("delta %g\n", delta);
  102. }
  103. else if ((t<=0.25)&&(x>-t+0.25)){delta = fabs(v_n_m - 1);}
  104. else if ((t>0.25)&&(fabs(x +(0.5*(t-0.25)))<0.00000001)) {delta = fabs(v_n_m -(4.0*x/(1-4.0*t)) );}
  105. else if ((t>0.25)&&(x<-0.5*(t-0.25))){delta = fabs(v_n_m);}
  106.  
  107. else if ((t>0.25)&&(x>-0.5*(t-0.25))) {delta = fabs (v_n_m - 1); };
  108. return delta;
  109. }
  110.  
  111.  
  112. void norma_c(double delta, double *max){
  113. if (fabs(delta) > *max) *max = fabs(delta);
  114. }
  115.  
  116.  
  117. void norma_l(double v_n_m, double * norma){
  118. *norma += fabs(v_n_m);
  119. }
  120.  
  121.  
  122. void print_matrix(double * matrix, int length){
  123. printf("\nMATRIX PRINT\n");
  124. int i,j;
  125. for (i = 0; i < length; i++)
  126. for (j = 0; j < length; j++){
  127. printf("%-4.5g ", *(matrix + i*length + j));
  128. }
  129. printf("\n");
  130. }
  131.  
  132.  
  133. void print_block_to_file(double * block, int length, FILE * fout){
  134. int i;
  135. for (i = 0; i < length; ++i)
  136. {
  137. fprintf(fout, "%-8g ", block[i]);
  138. }
  139. fprintf(fout, "\n");
  140. }
  141.  
  142. void begin_value(double * block, int length){
  143. int i;
  144. /*for (i = 0; i < length; ++i)
  145. {
  146. if (x_by_stepx(i) <= 0)
  147. block[i] = 1.0;
  148. else
  149. block[i] = 0.0;
  150. }*/
  151. for (i = 0; i < length; ++i)
  152. {
  153.  
  154.  
  155. if (x_by_stepx(i) <= 0)
  156. block[i] = 0.0;
  157. else if((x_by_stepx(i)>0)&&(x_by_stepx(i)<=0.25))
  158. {block[i] = 4.0*x_by_stepx(i);}
  159. else {block[i] = 1.0;}
  160. //printf("BLOCK[%g] = %g\n", x_by_stepx(i), block[i]);
  161. }
  162.  
  163.  
  164. }
  165.  
  166. void characteristic(double * block){
  167. /*block[0] = 1.0;
  168. block[m-1] = 0.0;*/
  169.  
  170. block[0] = 0.0;
  171. block[m-1] = 1.0;
  172. }
  173.  
  174. void characteristic_for_b(double * b){
  175. /*b[0] = 1-tao/(4*h);
  176. b[m-3] = 0.0;*/
  177. b[0] = 0.0;
  178. b[m-3] = 1+tao/(4*h);
  179. }
  180.  
  181. void next_block(double * block, int length, double step_by_t,//явный случай
  182. double * norma_DELTA_c_line, double * norma_delta_c_line, //[c]
  183. double * norma_DELTA_l_line, double * norma_delta_l_line //[l]
  184. ){
  185. int i;
  186. double x = x_begin;
  187. double del;
  188. *norma_DELTA_c_line = 0; *norma_delta_c_line = 0; //начальное максимальное значение [c]
  189. *norma_DELTA_l_line = 0; *norma_delta_l_line = 0; //начальное значение [l]
  190.  
  191. characteristic(block);
  192. for (i = 0; i < m-1; ++i)
  193. {
  194. block[i] = V_nplus1_m(block[i+1], block[i]);
  195. //block[i] = proverka(block[i]);
  196. del = delta(block[i], x, (step_by_t+1)*tao); // разность между численным решением и точным решением
  197. norma_c(del, norma_DELTA_c_line); //[c]
  198. norma_c(block[i], norma_delta_c_line); //[c]
  199. norma_l(del, norma_DELTA_l_line); //[l]
  200. norma_l(block[i], norma_delta_l_line);
  201. // printf("%g %g \n", del, x);
  202. x+=h;
  203. //if (del >= 1) { printf("t: %g\n ", tao*(step_by_t+1));};
  204. } //printf("----------------- \n");
  205.  
  206. // printf("NORMA 2 %g\n", *norma_delta_c_line);
  207.  
  208. }
  209.  
  210. void next_block_nelin(double * block, int length, double step_by_t,//явный случай
  211. double * norma_DELTA_c_line, double * norma_delta_c_line, //[c]
  212. double * norma_DELTA_l_line, double * norma_delta_l_line //[l]
  213. ){
  214. int i;
  215. double x = x_begin;
  216. double del;
  217. *norma_DELTA_c_line = 0; *norma_delta_c_line = 0; //начальное максимальное значение [c]
  218. *norma_DELTA_l_line = 0; *norma_delta_l_line = 0; //начальное значение [l]
  219.  
  220. characteristic(block);
  221. for (i = 0; i<m-1 ; ++i)
  222. {
  223. block[i] = V_nplus1_m_nelin(block[i+1], block[i]);
  224. //block[i] = proverka(block[i]);
  225. del = delta_nel(block[i], x, (step_by_t+1)*tao); // разность между численным решением и точным решением
  226. norma_c(del, norma_DELTA_c_line); //[c]
  227. norma_c(block[i], norma_delta_c_line); //[c]
  228. norma_l(del, norma_DELTA_l_line); //[l]
  229. norma_l(block[i], norma_delta_l_line);
  230. printf("%.15lf %g %lf \n", block[i], x, del);
  231. //if (del >= 1) { printf("t: %g\n ", tao*(step_by_t+1));};
  232. x+=h;
  233. //x = proverka(x);
  234. }
  235. printf("------------------");
  236. // printf("NORMA 2 %g\n", *norma_delta_c_line);
  237.  
  238. }
  239.  
  240. void matrix_b(double * b){
  241. int i;
  242. characteristic_for_b(b);
  243. //printf("b[%d]: %f\n", 0,b[0]);
  244. for (i = 1; i < m-1; i++)
  245. { if (x_by_stepx(i) <= 0)
  246. b[i] = 0.0;
  247. else if((x_by_stepx(i)>0)&&(x_by_stepx(i)<=0.25))
  248. {b[i] = 4*x_by_stepx(i);}
  249. else { b[i] = 1.0;}
  250. /*if (x_by_stepx(i) <= 0)
  251. b[i] = 1.0;
  252. else
  253. b[i] = 0.0;*/
  254. //printf("b[%d]: %f\n", i,b[i]);
  255. }
  256. //printf("b[%d]: %f\n", m-1,b[m-1]);
  257.  
  258.  
  259.  
  260.  
  261. }
  262.  
  263.  
  264.  
  265. void matrix_b1(double * b){
  266. int i;
  267. characteristic_for_b(b);
  268. for (i = 1; i < m-2; i++)
  269. {
  270. if ((x_by_stepx(i+2)==0)&&(x_by_stepx(i+1)!=0)){ // ОБЯЗАТЕЛЬНОЕ ПОПАДАНИЕ В 0
  271. b[i] = tao/8*(4*x_by_stepx(i+2));
  272. }
  273.  
  274. else if ((x_by_stepx(i+1)==0)&&(x_by_stepx(i)!=0)){
  275. b[i] = tao/8*(4*x_by_stepx(i+2)-16*x_by_stepx(i+1));
  276. }
  277.  
  278. else if ((x_by_stepx(i)==0)&&(x_by_stepx(i-1)!=0)){
  279. b[i] = 4*x_by_stepx(i)+tao/8*(4*x_by_stepx(i+2)-16*x_by_stepx(i+1)+24*x_by_stepx(i));
  280. }
  281.  
  282. else if ((x_by_stepx(i-1)==0)&&(x_by_stepx(i-2)!=0)){
  283. b[i] = 4*x_by_stepx(i+1)+tao/8*(4*x_by_stepx(i+2)-16*x_by_stepx(i+1)+24*x_by_stepx(i)-16*x_by_stepx(i-1));
  284. }
  285.  
  286. else if ((x_by_stepx(i-2)>0)&&(x_by_stepx(i+2)<=0.25)){
  287. b[i] = 4*x_by_stepx(i);
  288. }
  289.  
  290. else if ((x_by_stepx(i+2)>0.25)&&(x_by_stepx(i+1)<=0.25)){ // ОБЯЗАТЕЛЬНОЕ ПОПАДАНИЕ В 0
  291. b[i] = 4*x_by_stepx(i)+tao/8*(1-16*x_by_stepx(i+1)+24*x_by_stepx(i)-16*x_by_stepx(i-1)+4*x_by_stepx(i-2));
  292. }
  293.  
  294. else if ((x_by_stepx(i+1)>0.25)&&(x_by_stepx(i)<=0.25)){
  295. b[i] = 4*x_by_stepx(i)+tao/8*(-3+24*x_by_stepx(i)-16*x_by_stepx(i-1)+4*x_by_stepx(i-2));
  296. }
  297.  
  298. else if ((x_by_stepx(i)>0.25)&&(x_by_stepx(i-1)<=0.25)){
  299. b[i] = 1+tao/8*(3-16*x_by_stepx(i-1)+4*x_by_stepx(i-2));
  300. }
  301.  
  302. else if ((x_by_stepx(i-1)>0.25)&&(x_by_stepx(i-2)<=0.25)){
  303. b[i] = 1+tao/8*(-1+4*x_by_stepx(i-2));
  304. }
  305.  
  306. else if (x_by_stepx(i-2) > 0.25){
  307. b[i] = 1.0;
  308. }
  309. else
  310. b[i] = 0.0;
  311. }
  312. }
  313.  
  314.  
  315. void matrix_U( double U_1[m-2], double U_2[m-3], double b[m-2]){
  316. int i,j;
  317. double l=0; // начальное l
  318. for (i=0; i<m-2; i++){
  319. if (i==0) { //вырожденный случай
  320. U_1[i] = 1;
  321. U_2[i] = -tao/(4*h);
  322. b[i] = b[i];
  323. l = (tao/(4*h))/U_1[i];
  324. } else if (i!=m-3) {
  325. U_1[i] = 1 - l*U_2[i-1];
  326. U_2[i] = -tao/(4*h);
  327. b[i] = b[i] - l*b[i-1];
  328. l = (tao/(4*h))/U_1[i];
  329. } else {
  330. U_1[i] = 1 - l*U_2[i-1];
  331. b[i] = b[i] - l*b[i-1];
  332. }
  333.  
  334. }
  335. }
  336.  
  337. void matrix_V_hat(double U_1[m-2], double U_2[m-3], double b[m-2], double V_hat[m], double step_by_t,//неявный случай
  338. double * norma_DELTA_c_line, double * norma_delta_c_line, //[c]
  339. double * norma_DELTA_l_line, double * norma_delta_l_line){
  340. int i;
  341. double del;
  342. *norma_DELTA_c_line = 0; *norma_delta_c_line = 0; //начальное максимальное значение [c]
  343. *norma_DELTA_l_line = 0; *norma_delta_l_line = 0;
  344.  
  345.  
  346. characteristic(V_hat);
  347. //printf("V_hat[%d]: %g\n", m-1,V_hat[m-1]);
  348. for(i=m-2; i>=1; i--){
  349. if (i==m-2) V_hat[i] = b[i-1]/U_1[i-1];
  350. else {
  351. V_hat[i] = (b[i-1] - U_2[i-1]*V_hat[i+1])/U_1[i-1];
  352. V_hat[i] = proverka(V_hat[i]);
  353. }
  354.  
  355. del = delta(V_hat[i], x_by_stepx(i), (step_by_t+1)*tao); // разность между численным решением и точным решением
  356. //printf("del: %g\n", del);
  357. norma_c(del, norma_DELTA_c_line); //[c]
  358. norma_c(V_hat[i], norma_delta_c_line); //[c]
  359. norma_l(del, norma_DELTA_l_line); //[l]
  360. norma_l(V_hat[i], norma_delta_l_line);
  361. //printf("NORMA_DELTA_c_line %g\n", *norma_DELTA_c_line);
  362. //printf("%.10f %g\n",del, x_by_stepx(i));
  363. } //printf("--------------- \n");
  364.  
  365. }
  366.  
  367. void equal(double b[m-2], double V_hat[m]){
  368. int i;
  369. b[m-3] = V_hat[m-2]+tao/(4*h);
  370. b[m-3] = proverka(b[m-3]);
  371. for(i=0; i<m-3; i++){
  372. b[i] = V_hat[i+1];
  373. }
  374. }
  375.  
  376.  
  377.  
  378. void equal1(double b[m-2], double V_hat[m]){
  379. int i;
  380. b[0] = V_hat[1] +(tao/8.)*(V_hat[3] - 4*V_hat[2] + 6*V_hat[1] ) ;
  381. b[1] = V_hat[2] + (tao/8.)*(V_hat[4] - 4*V_hat[3] + 6*V_hat[2] - 4.0*V_hat[1] );
  382. for(i=2; i<m-4; i++){
  383. 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]);
  384. }
  385. 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]+1);
  386. b[m-3] = V_hat[m-2] + (tao/8.)*(V_hat[m-4] - 4*V_hat[m-3] + 6*V_hat[m-2]-3)+ tao/(4*h);
  387. }
  388.  
  389.  
  390. void main_loop_lin_simple(){
  391.  
  392. int i;
  393. double block[m];
  394. //по всей сетке [c]
  395. double norma_DELTA_c_line; double norma_delta_c_line; //по линиям
  396. //по всей сетке [l]
  397. double norma_DELTA_l_line; double norma_delta_l_line; //по линиям
  398. //FILE * fout = fopen("output_1.txt", "w+");
  399.  
  400. int length = (int)(sizeof(block)/sizeof(double));
  401. begin_value(block, length);
  402.  
  403. for (i = 0; i < n; ++i){
  404. next_block(block, length, i, &norma_DELTA_c_line, &norma_delta_c_line, &norma_DELTA_l_line, &norma_delta_l_line);
  405. /*norma_c(norma_DELTA_c_line, &norma_DELTA_c);
  406. norma_c(norma_delta_c_line, &norma_delta_c);
  407. norma_l(norma_DELTA_l_line, &norma_DELTA_l);
  408. norma_l(norma_delta_l_line, &norma_delta_l);*/
  409.  
  410. // print_block_to_file(block, length, fout);
  411.  
  412. }
  413.  
  414. printf("NORMA_DELTA_c %g\n", norma_DELTA_c_line);
  415. printf("NORMA_delta_c %g\n", norma_DELTA_c_line/norma_delta_c_line);
  416. printf("NORMA_DELTA_l %g\n", h*norma_DELTA_l_line);
  417. printf("NORMA_delta_l %g\n", norma_DELTA_l_line/norma_delta_l_line);
  418.  
  419. printf("tao h DeltaC DeltaL deltaC delta l\n");
  420. printf("%g & %g & %e & %e & %e & %e ", 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);
  421.  
  422. }
  423.  
  424. void main_loop_lin_nelin(){
  425. // double previous_block[m];
  426. int i;
  427. double block[m];
  428. //по всей сетке [c]
  429. double norma_DELTA_c_line; double norma_delta_c_line; //по линиям
  430. //по всей сетке [l]
  431. double norma_DELTA_l_line; double norma_delta_l_line; //по линиям
  432. FILE * fout = fopen("output_2.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_nelin(block, length, i, &norma_DELTA_c_line, &norma_delta_c_line, &norma_DELTA_l_line, &norma_delta_l_line);
  439.  
  440.  
  441. print_block_to_file(block, length, fout);
  442.  
  443. }
  444.  
  445. printf("NORMA_DELTA_c %g\n", norma_DELTA_c_line);
  446. printf("NORMA_delta_c %g\n", norma_DELTA_c_line/norma_delta_c_line);
  447. printf("NORMA_DELTA_l %g\n", h*norma_DELTA_l_line);
  448. printf("NORMA_delta_l %g\n", norma_DELTA_l_line/norma_delta_l_line);
  449.  
  450. printf("tao h DeltaC DeltaL deltaC delta l\n");
  451. printf("%g & %g & %e & %e & %e & %e ", 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);
  452.  
  453. }
  454.  
  455. //########################НЕЛИНЕЙНЫЙ СЛУЧАЙ###############################
  456. void begin_value_pristrel(double * block, int length){
  457. int i;
  458. for (i = 0; i < length; ++i){
  459. block[i] = 1.0;
  460. }
  461. }
  462.  
  463.  
  464. void matrix_J( double U_1[m-2], double U_2[m-3], double V_pr[m], double G_2[m-2]){
  465. int i,j;
  466. double l = 0; // íà÷àëüíîå l
  467. double a = 2*tao/(4*h);
  468.  
  469. for (i=0; i<m-2; i++){
  470. if (i==0) { //âûðîæäåííûé ñëó÷àé
  471. U_1[i] = 1;
  472. U_2[i] = -a*V_pr[2];
  473. G_2[i] = G_2[i];
  474. l = a*V_pr[0]/U_1[i]; // ????????
  475. } else if (i!=m-3) {
  476. U_1[i] = 1 - l*U_2[i-1];
  477. U_2[i] = -a*V_pr[i+2];
  478. G_2[i] = G_2[i] - l*G_2[i-1];
  479. l = a*V_pr[i+1]/U_1[i]; // ??????????
  480. } else {
  481. U_1[i] = 1 - l*U_2[i-1];
  482. G_2[i] = G_2[i] - l*G_2[i-1];
  483. }
  484. }
  485.  
  486.  
  487.  
  488. for (i = 0; i < m-3; ++i){
  489. //printf("U_1[%-2d]: %-13g, U_2[%-2d]: %-13g \n", i, U_1[i], i, U_2[i]);
  490. }
  491. }
  492.  
  493. void matrix_G(double V_pr[m], double G[m-2], double V[m]){
  494. int i;
  495. double a = 2*tao/(4*h);
  496.  
  497. for (i=0; i<m-2; i++){
  498. if (i==0) {
  499. G[0] = -(V_pr[1] - V[1] - a*(V_pr[2]*V_pr[2] - 0));
  500. } else if (i!=m-3) {
  501. G[i] = -(V_pr[i+1] - V[i+1] - a*(V_pr[i+2]*V_pr[i+2] - V_pr[i]*V_pr[i]));
  502. } else {
  503. G[i-3] = -(V_pr[i-2] - V[i-2] - a*(1 - V_pr[i-3]*V_pr[i-3]));
  504. }
  505.  
  506. }
  507. }
  508.  
  509. void matrix_G1(double V_pr[m], double G[m-2], double V[m]){
  510. int i;
  511. double a = 2*tao/(4*h);
  512.  
  513. for (i=0; i<m-2; i++){
  514. if (i==0) {
  515. G[0] = -(V_pr[1] - V[1] - a*(V_pr[2]*V_pr[2] - 0) -tao*0.125*(V[i+3]-4*V[i+2]+6*V[i]));
  516. } else if (i == 1)
  517. {
  518. 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]-4*V[i-1]));
  519. }
  520. else if ((i!=m-3)&&(i!=m-4) ){
  521. 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]-4*V[i-1]+V[i-2]));
  522. }
  523. else if (i== m-4)
  524. {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*(1-4*V[i+2]+6*V[i]-4*V[i-1]+V[i-2]));
  525. }
  526.  
  527. else { //i = m-3
  528. G[i-3] = -(V_pr[i-2] - V[i-2] - a*(1 - V_pr[i-3]*V_pr[i-3])- tao*0.125*(1-4+6*V[i]-4*V[i-1]+V[i-2]));
  529. }
  530.  
  531. }
  532. }
  533.  
  534.  
  535. void new_vector_X(double U_1[m-2], double U_2[m-3], double b[m-2], double X[m-2]){
  536.  
  537. int i;
  538. double del;
  539.  
  540. for(i=m-3; i>=0; i--){
  541. if (i == m-3) {
  542. X[i] = b[i]/U_1[i];
  543. }
  544. else {
  545. X[i] = (b[i] - U_2[i]*X[i+1])/U_1[i];
  546. }
  547. }
  548. }
  549.  
  550. void new_V_pr(double X[m-2], double V_pr[m]){
  551. int i;
  552. for (i = 0; i < m-2; ++i){
  553. V_pr[i+1] = V_pr[i+1] + X[i];
  554. }
  555. }
  556.  
  557. void equal_G2_G1(double G_1[m-2], double G_2[m-2]){
  558. int i;
  559. for (i = 0; i < m-2; ++i)
  560. {
  561. G_2[i] = G_1[i];
  562. }
  563. }
  564.  
  565. void equal_V_Vpr(double V[m], double V_pr[m]){
  566. int i;
  567. for (i = 0; i < m-2; ++i)
  568. {
  569. V[i] = V_pr[i];
  570. }
  571. }
  572.  
  573. double norma_G(double G[m-2]){
  574. double norma = 0;
  575. int i;
  576. for (i = 0; i < m-2; ++i){
  577. norma += fabs(G[i]);
  578. }
  579.  
  580. return norma*h;
  581. }
  582.  
  583.  
  584. int main(){
  585. printf("Size: %d x %d\n", m,n);
  586.  
  587. // main_loop_lin_simple();
  588. //main_loop_lin_nelin();
  589. //return 0;
  590. int i;
  591. //double A[m][m];
  592. double U_1[m-2]; double U_2[m-3]; double b[m-2]; double V_hat[m];
  593. 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; //чтобы зашло
  594.  
  595. int length = (int)(sizeof(V_hat)/sizeof(double));
  596.  
  597. //по всей сетке [c]
  598. double norma_DELTA_c_line; double norma_delta_c_line; //по линиям
  599. //по всей сетке [l]
  600. double norma_DELTA_l_line; double norma_delta_l_line; //по линиям
  601. FILE * fout = fopen("output_2.txt", "w+");
  602. FILE * fout_once = fopen("output_once.txt", "w");
  603. //##########################################
  604.  
  605.  
  606. i = 0;
  607.  
  608.  
  609. begin_value(V, length); // V === begin_value
  610.  
  611. for (i = 0; i < n; ++i)
  612. {
  613.  
  614. begin_value_pristrel(V_pr, length); // V_pr === begin_value
  615. V_pr[m-1] = 1;
  616. V_pr[0]=0;
  617.  
  618. while (norma_for_G > 0.001) {
  619.  
  620. //printf("%g - NORMAAAAAAAAAAAAAAAAAAAAAAAAAA GGGGGGGGGGGGG\n\n", norma_for_G);
  621.  
  622. matrix_G1(&V_pr, &G_1, &V);
  623. //print_block(G_1, length);
  624.  
  625. equal_G2_G1(&G_1, &G_2);
  626.  
  627. matrix_J(&U_1, &U_2, &V_pr, &G_2);
  628. /*printf("\n\n");
  629. print_block(U_1, m-2);
  630. printf("\n\n");
  631. print_block(U_2, m-3);
  632. printf("\n\n");*/
  633.  
  634. new_vector_X(&U_1, &U_2, &G_2, &X);
  635.  
  636. /*print_block(X, m-2);
  637. printf("\n\n");*/
  638.  
  639. new_V_pr(&X, &V_pr);
  640. //printf("XXXXXXX\n\n");
  641. //if (t_by_stept(i+1)>0.45) print_block(X, m-2);
  642.  
  643. matrix_G1(&V_pr, &G_1, &V); // изменение G для нового вектора пристреленного
  644.  
  645. printf("V_pr\n\n");
  646. print_block(V_pr, m);
  647. //if (t_by_stept(i+1)>0.45) printf("norma_g: %g \n", norma_for_G);
  648. norma_for_G = norma_G(&G_1);
  649. }
  650. norma_for_G = 1;
  651. equal_V_Vpr(&V, &V_pr);
  652. printf("t %g\n", t_by_stept(i+1));
  653. //if(t_by_stept(i+1)>0.25) break;
  654. }
  655.  
  656. return 0;
  657. //##############################################################
  658.  
  659. i = 0;
  660. matrix_b1(&b);
  661. for(i=0; i<n; i++){
  662. //printf("\nTao: %f, Step %d\n\n", t_by_stept(i+1), i);
  663.  
  664. matrix_U( &U_1, &U_2, &b);
  665. 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);
  666. equal1(&b, &V_hat);
  667. //printf("NORMA_DELTA_c_line %g\n", norma_DELTA_c_line);
  668. print_block_to_file(V_hat, length, fout);
  669.  
  670. //printf("NORMA_DELTA_c %g\n", norma_DELTA_c);
  671.  
  672. }
  673. //printf("NORMA_DELTA_c %g\n", norma_DELTA_c);
  674. //printf("NORMA_delta_c %g\n", norma_DELTA_c/norma_delta_c);
  675. //printf("NORMA_DELTA_l %g\n", h*norma_DELTA_l);
  676. printf("tao h DeltaC DeltaL deltaC delta l\n");
  677. printf("%g & %g & %.15lf & %e & %e & %e ", 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);
  678.  
  679. return 0;
  680. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!