API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Dec 14th, 2017
64
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 21.36 KB | None | 0 0
raw download clone embed print report
  1. #include<stdio.h>
  2. #include<stdlib.h>
  3. #include<string.h>
  4. #include<math.h>
  5.  
  6. //#define eps 10e-11
  7. #define PI 3.14159265358970323846
  8.  
  9. #define b21 1.0/5.0
  10.  
  11. #define b31 3.0/40.0
  12. #define b32 9.0/40.0
  13.  
  14. #define b41 44.0/45.0
  15. #define b42 -56.0/15.0
  16. #define b43 32.0/9.0
  17.  
  18. #define b51 19372.0/6561.0
  19. #define b52 -25360.0/2187.
  20. #define b53 64448.0/6561.0
  21. #define b54 -212.0/729.0
  22.  
  23. #define b61 9017.0/3168.0
  24. #define b62 -355.0/33.0
  25. #define b63 46732.0/5247.0
  26. #define b64 49.0/176.0
  27. #define b65 -5103.0/18656.0
  28.  
  29. #define b71 35.0/384.0
  30. #define b72 0.
  31. #define b73 500.0/1113.0
  32. #define b74 125.0/192.0
  33. #define b75 -2187.0/6784.0
  34. #define b76 11.0/84.0
  35.  
  36. #define c1 35.0/384.0
  37. #define c2 0.
  38. #define c3 500.0/1113.0
  39. #define c4 125.0/192.0
  40. #define c5 -2187.0/6784.0
  41. #define c6 11.0/84.0
  42. #define c7 0.
  43.  
  44. #define d1 5179.0/57600.0
  45. #define d2 0.
  46. #define d3 7571.0/16695.0
  47. #define d4 393.0/640.0
  48. #define d5 -92097.0/339200.0
  49. #define d6 187.0/2100.0
  50. #define d7 1.0/40.0
  51.  
  52.  
  53. // Функции (5 штук)
  54. double fx(double x,double y,double z, double u){
  55.     return y;
  56. }
  57. double fy(double x,double y,double z,double u,double alpha){ //дано начальное знаечние
  58.     return u - x * exp( -alpha*x);
  59. }
  60. double fz(double x,double y,double z,double u,double alpha){ //дано начальное знаечние
  61.     return -u*exp(-alpha*x) + alpha*u*x*exp(-alpha*x);
  62. }
  63. double fu(double x,double y,double z,double u){
  64.     return z;
  65. }
  66. double fv(double x,double y,double z,double u){ // наша функция есть производная этого интеграла
  67.     return u*u;
  68. }
  69.  
  70.  
  71.  
  72. double Xnewvar(double h,double x,double kx1,double kx2,double kx3,double kx4,double kx5,double kx6,double kx7){
  73.     return x+(c1*kx1+c2*kx2+c3*kx3+c4*kx4+c5*kx5+c6*kx6+c7*kx7);
  74. }
  75. double Ynewvar(double h,double y,double ky1,double ky2,double ky3,double ky4,double ky5,double ky6,double ky7){
  76.     return y+(c1*ky1+c2*ky2+c3*ky3+c4*ky4+c5*ky5+c6*ky6+c7*ky7);
  77. }
  78. double Znewvar(double h,double z,double kz1,double kz2,double kz3,double kz4,double kz5,double kz6,double kz7){
  79.     return z+(c1*kz1+c2*kz2+c3*kz3+c4*kz4+c5*kz5+c6*kz6+c7*kz7);
  80. }
  81. double Unewvar(double h,double u,double ku1,double ku2,double ku3,double ku4,double ku5,double ku6,double ku7){
  82.     return u+(c1*ku1+c2*ku2+c3*ku3+c4*ku4+c5*ku5+c6*ku6+c7*ku7);
  83. }
  84. double Vnewvar(double h,double v,double kv1,double kv2,double kv3,double kv4,double kv5,double kv6,double kv7){
  85.     return v+(c1*kv1+c2*kv2+c3*kv3+c4*kv4+c5*kv5+c6*kv6+c7*kv7);
  86. }
  87.  
  88.  
  89. double Xnewvartilde(double h,double x,double kx1,double kx2,double kx3,double kx4,double kx5,double kx6,double kx7){
  90.     return x+(d1*kx1+d2*kx2+d3*kx3+d4*kx4+d5*kx5+d6*kx6+d7*kx7);
  91. }
  92. double Ynewvartilde(double h,double y,double ky1,double ky2,double ky3,double ky4,double ky5,double ky6,double ky7){
  93.     return y+(d1*ky1+d2*ky2+d3*ky3+d4*ky4+d5*ky5+d6*ky6+d7*ky7);
  94. }
  95. double Znewvartilde(double h,double z,double kz1,double kz2,double kz3,double kz4,double kz5,double kz6,double kz7){
  96.     return z+(d1*kz1+d2*kz2+d3*kz3+d4*kz4+d5*kz5+d6*kz6+d7*kz7);
  97. }
  98. double Unewvartilde(double h,double u,double ku1,double ku2,double ku3,double ku4,double ku5,double ku6,double ku7){
  99.     return u+(d1*ku1+d2*ku2+d3*ku3+d4*ku4+d5*ku5+d6*ku6+d7*ku7);
  100. }
  101. double Vnewvartilde(double h,double v,double kv1,double kv2,double kv3,double kv4,double kv5,double kv6,double kv7){
  102.     return v+(d1*kv1+d2*kv2+d3*kv3+d4*kv4+d5*kv5+d6*kv6+d7*kv7);
  103. }
  104.  
  105.  
  106. double delta(double x,double y,double z,double u, double v,double x1,double y1,double z1,double u1, double v1){
  107.     return sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1)+(u-u1)*(u-u1)+(v-v1)*(v-v1));
  108. }
  109.  
  110. // Максимальное, минимальное, норма,
  111. double max(double a,double b){
  112.     if (a>b) return a;
  113.     else return b;
  114. }
  115. double min(double a,double b){
  116.     if (a<b) return a;
  117.     else return b;
  118. }
  119. double mod(double a,double b){
  120.     return sqrt(a*a+b*b);
  121. }
  122.  
  123. // [на самом деле это Рунге Кутта 5 порядка]
  124. double err(double *t,double *x,double *y,double *z,double *u, double *v,
  125.            double *h,double facmax, FILE*out,
  126.            double eps,double T,double *k,double *dx,double *dy, double alpha)
  127. {
  128.     // facmax - хрень для выбора шага
  129.     double tol=eps;
  130.     double p = 1.0/8.0;
  131.     double kx1,kx2,kx3,kx4,kx5,kx6,kx7;
  132.     double ky1,ky2,ky3,ky4,ky5,ky6,ky7;
  133.     double kz1,kz2,kz3,kz4,kz5,kz6,kz7;
  134.     double ku1,ku2,ku3,ku4,ku5,ku6,ku7;
  135.     double kv1,kv2,kv3,kv4,kv5,kv6,kv7;
  136.  
  137.     double x1,y1,z1,u1,v1;
  138.     double xx,yy,zz,uu,vv;
  139.     double d,hnew;
  140.     double ddx, ddy;
  141.     //fprintf(out,"X%g %g %g %g\n",*x,*y,*z,*u);
  142.  
  143.     kx1=*h*fx(*x,*y,*z,*u);
  144.     ky1=*h*fy(*x,*y,*z,*u, alpha);
  145.     kz1=*h*fz(*x,*y,*z,*u, alpha);
  146.     ku1=*h*fu(*x,*y,*z,*u);
  147.     kv1=*h*fv(*x,*y,*z,*u);
  148.  
  149.     kx2=*h*fx(*x+b21*kx1,*y+b21*ky1,*z+b21*kz1,*u+b21*ku1);
  150.     ky2=*h*fy(*x+b21*kx1,*y+b21*ky1,*z+b21*kz1,*u+b21*ku1, alpha);
  151.     kz2=*h*fz(*x+b21*kx1,*y+b21*ky1,*z+b21*kz1,*u+b21*ku1, alpha);
  152.     ku2=*h*fu(*x+b21*kx1,*y+b21*ky1,*z+b21*kz1,*u+b21*ku1);
  153.     kv2=*h*fv(*x+b21*kx1,*y+b21*ky1,*z+b21*kz1,*u+b21*ku1);
  154.  
  155.  
  156.     kx3=*h*fx(*x+b31*kx1+b32*kx2,*y+b31*ky1+b32*ky2,*z+b31*kz1+b32*kz2,*u+b31*ku1+b32*ku2);
  157.     ky3=*h*fy(*x+b31*kx1+b32*kx2,*y+b31*ky1+b32*ky2,*z+b31*kz1+b32*kz2,*u+b31*ku1+b32*ku2, alpha);
  158.     kz3=*h*fz(*x+b31*kx1+b32*kx2,*y+b31*ky1+b32*ky2,*z+b31*kz1+b32*kz2,*u+b31*ku1+b32*ku2, alpha);
  159.     ku3=*h*fu(*x+b31*kx1+b32*kx2,*y+b31*ky1+b32*ky2,*z+b31*kz1+b32*kz2,*u+b31*ku1+b32*ku2);
  160.     kv3=*h*fv(*x+b31*kx1+b32*kx2,*y+b31*ky1+b32*ky2,*z+b31*kz1+b32*kz2,*u+b31*ku1+b32*ku2);
  161.  
  162.  
  163.     kx4=*h*fx(*x+b41*kx1+b42*kx2+b43*kx3,
  164.             *y+b41*ky1+b42*ky2+b43*ky3,
  165.             *z+b41*kz1+b42*kz2+b43*kz3,
  166.             *u+b41*ku1+b42*ku2+b43*ku3);
  167.     ky4=*h*fy(*x+b41*kx1+b42*kx2+b43*kx3,
  168.             *y+b41*ky1+b42*ky2+b43*ky3,
  169.             *z+b41*kz1+b42*kz2+b43*kz3,
  170.             *u+b41*ku1+b42*ku2+b43*ku3, alpha);
  171.     kz4=*h*fz(*x+b41*kx1+b42*kx2+b43*kx3,
  172.             *y+b41*ky1+b42*ky2+b43*ky3,
  173.             *z+b41*kz1+b42*kz2+b43*kz3,
  174.             *u+b41*ku1+b42*ku2+b43*ku3, alpha);
  175.     ku4=*h*fu(*x+b41*kx1+b42*kx2+b43*kx3,
  176.             *y+b41*ky1+b42*ky2+b43*ky3,
  177.             *z+b41*kz1+b42*kz2+b43*kz3,
  178.             *u+b41*ku1+b42*ku2+b43*ku3);
  179.     kv4=*h*fv(*x+b41*kx1+b42*kx2+b43*kx3,
  180.             *y+b41*ky1+b42*ky2+b43*ky3,
  181.             *z+b41*kz1+b42*kz2+b43*kz3,
  182.             *u+b41*ku1+b42*ku2+b43*ku3);
  183.  
  184.     kx5=*h*fx(*x+b51*kx1+b52*kx2+b53*kx3+b54*kx4,
  185.             *y+b51*ky1+b52*ky2+b53*ky3+b54*ky4,
  186.             *z+b51*kz1+b52*kz2+b53*kz3+b54*kz4,
  187.             *u+b51*ku1+b52*ku2+b53*ku3+b54*ku4);
  188.     ky5=*h*fy(*x+b51*kx1+b52*kx2+b53*kx3+b54*kx4,
  189.             *y+b51*ky1+b52*ky2+b53*ky3+b54*ky4,
  190.             *z+b51*kz1+b52*kz2+b53*kz3+b54*kz4,
  191.             *u+b51*ku1+b52*ku2+b53*ku3+b54*ku4, alpha);
  192.     kz5=*h*fz(*x+b51*kx1+b52*kx2+b53*kx3+b54*kx4,
  193.             *y+b51*ky1+b52*ky2+b53*ky3+b54*ky4,
  194.             *z+b51*kz1+b52*kz2+b53*kz3+b54*kz4,
  195.             *u+b51*ku1+b52*ku2+b53*ku3+b54*ku4, alpha);
  196.     ku5=*h*fu(*x+b51*kx1+b52*kx2+b53*kx3+b54*kx4,
  197.             *y+b51*ky1+b52*ky2+b53*ky3+b54*ky4,
  198.             *z+b51*kz1+b52*kz2+b53*kz3+b54*kz4,
  199.             *u+b51*ku1+b52*ku2+b53*ku3+b54*ku4);
  200.     kv5=*h*fv(*x+b51*kx1+b52*kx2+b53*kx3+b54*kx4,
  201.             *y+b51*ky1+b52*ky2+b53*ky3+b54*ky4,
  202.             *z+b51*kz1+b52*kz2+b53*kz3+b54*kz4,
  203.             *u+b51*ku1+b52*ku2+b53*ku3+b54*ku4);
  204.  
  205.     kx6=*h*fx(*x+b61*kx1+b62*kx2+b63*kx3+b64*kx4+b65*kx5,
  206.             *y+b61*ky1+b62*ky2+b63*ky3+b64*ky4+b65*ky5,
  207.             *z+b61*kz1+b62*kz2+b63*kz3+b64*kz4+b65*kz5,
  208.             *u+b61*ku1+b62*ku2+b63*ku3+b64*ku4+b65*ku5);
  209.     ky6=*h*fy(*x+b61*kx1+b62*kx2+b63*kx3+b64*kx4+b65*kx5,
  210.             *y+b61*ky1+b62*ky2+b63*ky3+b64*ky4+b65*ky5,
  211.             *z+b61*kz1+b62*kz2+b63*kz3+b64*kz4+b65*kz5,
  212.             *u+b61*ku1+b62*ku2+b63*ku3+b64*ku4+b65*ku5, alpha);
  213.     kz6=*h*fz(*x+b61*kx1+b62*kx2+b63*kx3+b64*kx4+b65*kx5,
  214.             *y+b61*ky1+b62*ky2+b63*ky3+b64*ky4+b65*ky5,
  215.             *z+b61*kz1+b62*kz2+b63*kz3+b64*kz4+b65*kz5,
  216.             *u+b61*ku1+b62*ku2+b63*ku3+b64*ku4+b65*ku5, alpha);
  217.     ku6=*h*fu(*x+b61*kx1+b62*kx2+b63*kx3+b64*kx4+b65*kx5,
  218.             *y+b61*ky1+b62*ky2+b63*ky3+b64*ky4+b65*ky5,
  219.             *z+b61*kz1+b62*kz2+b63*kz3+b64*kz4+b65*kz5,
  220.             *u+b61*ku1+b62*ku2+b63*ku3+b64*ku4+b65*ku5);
  221.     kv6=*h*fv(*x+b61*kx1+b62*kx2+b63*kx3+b64*kx4+b65*kx5,
  222.             *y+b61*ky1+b62*ky2+b63*ky3+b64*ky4+b65*ky5,
  223.             *z+b61*kz1+b62*kz2+b63*kz3+b64*kz4+b65*kz5,
  224.             *u+b61*ku1+b62*ku2+b63*ku3+b64*ku4+b65*ku5);
  225.  
  226.     kx7=*h*fx(*x+b71*kx1+b72*kx2+b73*kx3+b74*kx4+b75*kx5+b76*kx6,
  227.             *y+b71*ky1+b72*ky2+b73*ky3+b74*ky4+b75*ky5+b76*ky6,
  228.             *z+b71*kz1+b72*kz2+b73*kz3+b74*kz4+b75*kz5+b76*kz6,
  229.             *u+b71*ku1+b72*ku2+b73*ku3+b74*ku4+b75*ku5+b76*ku6);
  230.     ky7=*h*fy(*x+b71*kx1+b72*kx2+b73*kx3+b74*kx4+b75*kx5+b76*kx6,
  231.             *y+b71*ky1+b72*ky2+b73*ky3+b74*ky4+b75*ky5+b76*ky6,
  232.             *z+b71*kz1+b72*kz2+b73*kz3+b74*kz4+b75*kz5+b76*kz6,
  233.             *u+b71*ku1+b72*ku2+b73*ku3+b74*ku4+b75*ku5+b76*ku6, alpha);
  234.     kz7=*h*fz(*x+b71*kx1+b72*kx2+b73*kx3+b74*kx4+b75*kx5+b76*kx6,
  235.             *y+b71*ky1+b72*ky2+b73*ky3+b74*ky4+b75*ky5+b76*ky6,
  236.             *z+b71*kz1+b72*kz2+b73*kz3+b74*kz4+b75*kz5+b76*kz6,
  237.             *u+b71*ku1+b72*ku2+b73*ku3+b74*ku4+b75*ku5+b76*ku6, alpha);
  238.     ku7=*h*fu(*x+b71*kx1+b72*kx2+b73*kx3+b74*kx4+b75*kx5+b76*kx6,
  239.             *y+b71*ky1+b72*ky2+b73*ky3+b74*ky4+b75*ky5+b76*ky6,
  240.             *z+b71*kz1+b72*kz2+b73*kz3+b74*kz4+b75*kz5+b76*kz6,
  241.             *u+b71*ku1+b72*ku2+b73*ku3+b74*ku4+b75*ku5+b76*ku6);
  242.     kv7=*h*fv(*x+b71*kx1+b72*kx2+b73*kx3+b74*kx4+b75*kx5+b76*kx6,
  243.             *y+b71*ky1+b72*ky2+b73*ky3+b74*ky4+b75*ky5+b76*ky6,
  244.             *z+b71*kz1+b72*kz2+b73*kz3+b74*kz4+b75*kz5+b76*kz6,
  245.             *u+b71*ku1+b72*ku2+b73*ku3+b74*ku4+b75*ku5+b76*ku6);
  246.  
  247.     xx=Xnewvar(*h,*x,kx1,kx2,kx3,kx4,kx5,kx6,kx7);
  248.     yy=Ynewvar(*h,*y,ky1,ky2,ky3,ky4,ky5,ky6,ky7);
  249.     zz=Znewvar(*h,*z,kz1,kz2,kz3,kz4,kz5,kz6,kz7);
  250.     uu=Unewvar(*h,*u,ku1,ku2,ku3,ku4,ku5,ku6,ku7);
  251.     vv=Vnewvar(*h,*v,kv1,kv2,kv3,kv4,kv5,kv6,kv7);
  252.  
  253.  
  254.     x1=Xnewvartilde(*h,*x,kx1,kx2,kx3,kx4,kx5,kx6,kx7);
  255.     y1=Ynewvartilde(*h,*y,ky1,ky2,ky3,ky4,ky5,ky6,ky7);
  256.     z1=Znewvartilde(*h,*z,kz1,kz2,kz3,kz4,kz5,kz6,kz7);
  257.     u1=Unewvartilde(*h,*u,ku1,ku2,ku3,ku4,ku5,ku6,ku7);
  258.     v1=Vnewvartilde(*h,*v,kv1,kv2,kv3,kv4,kv5,kv6,kv7);
  259.  
  260.     d = delta(xx,yy,zz,uu,vv,x1,y1,z1,u1,v1);
  261.  
  262.     hnew=(*h)*min(facmax,max(0.3,pow(tol/d,p))); // ***MAGIC*** ???? YES! tol = epsilon | d - ошибка
  263.  
  264.     if (d < eps) {d = eps;}
  265.  
  266.     *t+=*h;
  267.     *x=xx;  *y=yy; *z=zz; *u=uu;
  268.     //printf("ok2 t: %e x: %e y: %e z: %e u: %e\n",*t,kx1,ky1,kz1,ku1); fprintf(out,"t: %e x: %e y: %e z: %e u: %e\n",*t,*x,*y,*z,*u);
  269.     //printf("X%g %g %g %g\n",*x,*y,*z,*u);
  270.  
  271.     if (hnew>10e-5) hnew=10e-5; // Шаг не больше чем 10e-5
  272.     if (d > eps){return err(&*t,&*x,&*y,&*z,&*u,&*v,&*h,1.0,out,eps,T,&*k,&*dx,&*dy,alpha);} // вызывваем себя [facmax=1.0]
  273.     else
  274.     if ((*t+*h)>T) {*h=T-*t; return err(&*t,&*x,&*y,&*z,&*u,&*v,&*h,facmax,out,eps,T,&*k,&*dx,&*dy,alpha);}
  275.     else
  276.         if ((*t)==T){
  277.             //printf("%f %f %f %f %f %f\n",*t,*x,*y,*z,*u,*v);
  278.             //fprintf(out,"%f %f %f %f %f %f\n",*t,*x,*y,*z,*u,*v);
  279.             ++*k;
  280.             return d;
  281.             //printf("%f %f %f %f %f\n",*h,*u,*y,*x,(*h)*((*u)*(*u)-(*y)*(*y)-(*x)*(*x)));
  282.             //return (*h)*((*u)*(*u)-(*y)*(*y)-(*x)*(*x));
  283.         }
  284.         else{
  285.             *h=hnew;
  286.             //printf("%f %f %f %f %f %f\n",*t,*x,*y,*z,*u,*v);
  287.             //fprintf(out,"%f %f %f %f %f %f\n",*t,*x,*y,*z,*u,*v);
  288.             ddx=fabs(sin(*t)-*x); // не интерисует
  289.             ddy=fabs(cos(*t)-*y); // не интерисует
  290.             if (ddx>*dx) *dx=ddx; // не интерисует
  291.             if (ddy>*dy) *dy=ddy; // не интерисует
  292.             ++*k;
  293.             return  d;
  294.             // printf("%f %f %f %f %f\n",*h,*u,*y,*x,(*h)*((*u)*(*u)-(*y)*(*y)-(*x)*(*x)));
  295.             // return (*h)*((*u)*(*u)-(*y)*(*y)-(*x)*(*x));
  296.         }
  297. }
  298.  
  299. // [Рунге кутта]
  300. double runge(double x0,double u0,double *X,double *U, FILE*out)
  301. {
  302.     double t=0;
  303.     double h;
  304.     double x,y,z,u,v;
  305.     double eps;
  306.     double facmax=2.0; // ?
  307.     double er;
  308.     double del;
  309.     double T;
  310.     double I=0;
  311.         double alpha=25.0;
  312.  
  313.     double dx7,dx9,dx11;
  314.     double dy7,dy9,dy11;
  315.     double del7,del9,del11;
  316.     double x7,y7,z7,u7,v7;
  317.     double x9,y9,z9,u9,v9;
  318.     double x11,y11,z11,u11,v11;
  319.     double k,dx,dy;
  320.  
  321.     k=0;  // Количесвто шагов в начале
  322.     dx=0;
  323.     dy=0;
  324.  
  325.     T=PI*0.5;// t конечное
  326.  
  327.     // У меня было заданы начальные условия x и y для z и u используется метод стрельбы
  328.     eps=10e-15;
  329.     h=eps;
  330.     y=0; z=0; // Катины начальные условия
  331.     x=x0; u=u0; // Катины неизвестные нач. условия
  332.     v=0; //начальное значение интеграла?
  333.  
  334.     facmax=2.0; er=0; del=0; t=0; // обнулечние начальных
  335.     double d=0; // то что возвращает err
  336.  
  337.     while (t<T){
  338.         d = err(&t,&x,&y,&z,&u,&v,&h,facmax,out,eps,T,&k,&dx,&dy,alpha); // z и u - не знаем в начале, x и y - значем в начале
  339.         fprintf(out,"%.20lf %lf\n",t,u);
  340.         //I+=0.5*(u*u-y*y-x*x - (uprev*uprev - yprev*yprev - xprev*xprev))*(t-tprev);
  341.         /*Tm=(t-tprev)*0.5;
  342.         xm=xprev;
  343.         ym=yprev;
  344.         zm=zprev;
  345.         um=uprev;
  346.         tm=tprev;
  347.         d=err(&tm,&xm,&ym,&zm,&um,&Tm,facmax,out,eps,T,&k,&dx,&dy);
  348.         I+=0.5*(um*um - ym*ym - xm*xm)*(t-tprev);*/
  349.         }
  350.  
  351.     //А в конце были заданы значения x и u
  352.     *X=x;
  353.     *U=u;
  354.  
  355.     //printf("INTEGRAL: %.20lf\n", v);
  356.     return v;
  357. }
  358.  
  359.  
  360. // [Метод стрельбы]
  361. double sp(double *x0,double *u0,double *X,double *U, FILE*out) // x0 и u0 - не даны в начале
  362. {
  363.     double gamma=1; // Начальная гамма
  364.     double x1,u1;
  365.     double eps=10e-12;
  366.     double a,b,c,d;
  367.     double X11,X12,X21,X22; // нужно для подсчета матрицы
  368.     double X1,X2;
  369.     double x,u;
  370.     double xx,uu;
  371.     double det;
  372.  
  373.  
  374.     double XX=0.8; double UU=0; // известные конечные значения для x и u
  375.     //X1=X(T) X2=u(T)-1
  376.     fprintf(out,"00\n");
  377.     runge(*x0,*u0,&x,&u,out);
  378.     X1=x-XX; X2=u-UU; // Функции невязки
  379.     xx=x-XX; uu=u-UU; // Функции невязки
  380.     //fprintf(out,"eps  + 0\n");
  381.     runge(*x0+eps,*u0,&x,&u,out); // подсчет производной
  382.     X11=x-XX; X21=u-UU;
  383.     //fprintf(out,"0 + eps\n");
  384.     runge(*x0,*u0+eps,&x,&u,out); // подсчет производной
  385.     X12=x-XX; X22=u-UU;
  386.  
  387.     // элементы производной матрицы
  388.     a=(X11-X1)/eps;
  389.     b=(X12-X1)/eps;
  390.     c=(X21-X2)/eps;
  391.     d=(X22-X2)/eps;
  392.  
  393.     det=a*d-b*c; // опредилитель матрицы производной функции невязки
  394.  
  395.     //printf(out,"x0 %g %g\n",det,det);
  396.     if (det < eps && det > -eps) det = eps;
  397.     //fprintf(out,"x0 %g %g\n",det,det);
  398.     double h1=(d*X1-b*X2)/det;  
  399.     double h2=-(X1*c-X2*a)/det;  
  400.     x1=*x0-gamma*h1;
  401.     u1=*u0-gamma*h2;
  402.     printf("H1 = %f, H2 = %f\n", h1,h2);
  403.     while (mod(X1,X2) >= mod(xx,uu)){ //переходим к следующей гамме
  404.         gamma*=0.5;
  405.         if (gamma < eps) break;
  406.         x1=*x0-gamma*h1;
  407.         u1=*u0-gamma*h2;    
  408.         //fprintf(out,"chikl %f %f \n",x1,u1);
  409.             x=0; u=0;
  410.             runge(x1,u1,&x,&u,out);
  411.         printf("        X1 = %g U1 = %g\n", x1,u1);
  412.         X1=x-XX; X2=u-UU;
  413.     }
  414.  
  415.     // теперь гамма хорошая
  416.     *x0=x1;  // нашли хорошие начальные условия
  417.     *u0=u1;  // нашли хорошие начальные условия
  418.     //runge(*x0,*u0,&x,&u,out);
  419.  
  420.     X1=x-XX; X2=u-UU; // окончательная функция невязки?
  421.     *X=X1+XX;
  422.     *U=X2+UU;
  423.  
  424.     return 0;
  425. }
  426.  
  427. void pristrel(FILE*out)
  428. {
  429.  
  430.     double x,u;
  431.     double x0,u0;
  432.     double d=0;
  433.     double XX,UU; // в конце заданные x и u
  434.     XX=0.8; UU=0; // их значения
  435.     double v=0;   //
  436.     double eps=10e-8;
  437.  
  438.     //начальные параметры пристрелки
  439.     //x0=-193.9873696474; u0=-1246.4150443086;//T=20;
  440.     //x0=36.191649591711020; u0=78.592331658908350;//T=10;
  441.    
  442.     x0=0;  u0=10;   // начальные условия (которые будем менять)
  443.     x=0.8; u=0;  // конечные условия  
  444.  
  445.     //X1=X(T) X2=u(T)-1
  446.     runge(x0,u0,&x,&u,out);
  447.     d=mod(x-XX,u-UU);
  448.     printf("d: %g x: %g u: %g x0: %g u0: %g.20\n",d,x,u,x0,u0);
  449.     //fprintf(out,"oi%1.20f %1.20f %1.20f %1.20f\n",x,u,x0,u0);
  450.     printf("%f\n",d);
  451.     //printf("%f\n",d);
  452.     while(d > eps){
  453.         sp(&x0,&u0,&x,&u,out);
  454.         //runge(x0,u0,&x,&u,out);
  455.         d = mod(x-XX,u-UU);
  456.         //fprintf(out,"%f\n",d);
  457.         printf("d: %.20lf x: %g u: %g x0: %g u0: %lf \n",d,x,u,x0,u0);
  458.         /*fprintf(out,"%e %e %e %e\n",x,u,x0,u0);*/
  459.     }
  460.  
  461.     v = runge(x0,u0,&x,&u,out); // значение интеграла
  462.  
  463.     printf("x0: %1.15f u0: %1.15f\n",x0,u0);
  464.     printf("x: %1.15f u: %1.15f",x,u);
  465.  
  466.     /*fprintf(out,"Integral: %e\n",v);
  467.     fprintf(out,"x0: %e u0: %1e\n",x0,u0);
  468.     fprintf(out,"x: %e u: %e",x,u);*/
  469.  
  470. }
  471.  
  472. double zapuskrunge(FILE*out)
  473. {
  474.     double x,u;
  475.     double x0,u0;
  476.     double I;
  477.  
  478.     x0=-0.079004795578223; u0=0.954478985578936; //T=Pi/2;  // как я понимаю начальные усовия уже посчитали
  479.  
  480.     //z0=36.191649591711020; u0=78.592331658908350;//T=10;
  481.     //z0=-100; u0=-100;
  482.     x=0;  u=0;
  483.  
  484.     I=runge(x0,u0,&x,&u,out);
  485.  
  486.     printf("x0: %1.10f u0: %1.10f\n",x0,u0);
  487.     printf("x: %1.10f u: %1.10f",x,u);
  488.     printf("\n Integral %f\n",I);
  489.     printf("x0: %1.10f u0: %1.10f\n",x0,u0);
  490.     printf("x: %1.10f u: %1.10f",x,u);
  491.     //fprintf(out,"\n Integral %f\n",I);
  492.  
  493.     //fprintf(out,"\n Integral %f\n",I);
  494.  
  495.     //fprintf(out,"z0: %1.10f u0: %1.10f\n",z0,u0);
  496.     //fprintf(out,"x: %1.10f u: %1.10f",x,u);
  497.     //
  498.     return 0;
  499. }
  500.  
  501. double test(FILE*out) // для различных epsilon и выводим числа рунге
  502. {
  503.     double t=0;
  504.     double h;
  505.     double x,y,z,u,v;
  506.     double eps;
  507.     double facmax=2.0;
  508.     double er=0;
  509.     double del=0;
  510.     double T=PI;
  511.     double dx,dy;
  512.     double k; //int k;
  513.     double dx7,dx9,dx11;
  514.     double dy7,dy9,dy11;
  515.     int k7,k9,k11;
  516.     double del7,del9,del11;
  517.     double x7,y7;
  518.     double x9,y9;
  519.     double x11,y11;
  520.     int R=1;
  521.     double x0,y0,z0,u0;
  522.     x0=0;y0=1;z0=0;u0=0;
  523.     double gamma=1.;
  524.         double alpha=5;
  525.  
  526.     T=PI;
  527.     while (R!=10e6){
  528.  
  529.         T=R*PI;
  530.  
  531.         eps=10e-7; k=0; dx=0; dy=0;
  532.         h=eps;  x=x0; y=y0; z=z0; u=u0; facmax=2.0; er=0; del=0; t=0;
  533.         //printf("ok%e %e %e %e %e\n",x,y,z,u,h);
  534.         while (t<gamma*T){
  535.             del+=err(&t,&x,&y,&z,&u,&v,&h,facmax,out,eps,gamma*T,&k,&dx,&dy,alpha);
  536.             //printf("ok%e %e %e %e %e\n",x,y,z,u,h);
  537.         }
  538.         del7=del;k7=k;
  539.         dx7=dx;dy7=dy;
  540.         x7=x; y7=y;
  541.  
  542.  
  543.         eps=10e-9;  k=0; dx=0; dy=0;
  544.         h=eps;x=x0; y=y0; z=z0; u=u0; facmax=2.0; er=0; del=0; t=0;
  545.         while (t<gamma*T){
  546.             del+=err(&t,&x,&y,&z,&u,&v,&h,facmax,out,eps,gamma*T,&k,&dx,&dy,alpha);
  547.         }
  548.         del9=del;k9=k;
  549.         dx9=dx;dy9=dy;
  550.         x9=x;  y9=y;
  551.  
  552.         eps=10e-11;  k=0; dx=0; dy=0;
  553.         h=eps; x=x0; y=y0; z=z0; u=u0; facmax=2.0; er=0; del=0; t=0;
  554.         while (t<gamma*T){
  555.             del+=err(&t,&x,&y,&z,&u,&v,&h,facmax,out,eps,gamma*T,&k,&dx,&dy,alpha);
  556.         }
  557.         del11=del; k11=k;
  558.         dx11=dx; dy11=dy;
  559.         x11=x;  y11=y;
  560.  
  561.         //k-steps
  562.  
  563.         printf("T: %d Pi\n\n",R);
  564.         printf("eps 10e-7 k: %d, |x(T)|: %e  |y(T)-cos(T)|: %e dx:%e dy: %e  delta(T): %e\n\n",k7,fabs(x7),fabs(y7-cos(T)),dx7,dy7,del7);
  565.         printf("eps 10e-9 k: %d, |x(T)|: %e  |y(T)-cos(T)|: %e dx:%e dy: %e  delta(T): %e\n\n",k9,fabs(x9),fabs(y9-cos(T)),dx9,dy9,del9);
  566.         printf("eps 10e-11 k: %d, |x(T)|: %e  |y(T)-cos(T)|: %e dx:%e dy: %e  delta(T): %e\n\n",k11,fabs(x11),fabs(y11-cos(T)),dx11,dy11,del11);
  567.         printf("Runge number x for T: %g\n\n", fabs((x7-x9)/(x9-x11)));
  568.         printf("Runge number y for T: %g\n\n", fabs((y7-y9)/(y9-y11)));
  569.  
  570.         fprintf(out,"T: %d Pi\n",R);
  571.         fprintf(out,"eps 10e-7 k: %d, |x(T)|: %e  |y(T)-cos(T)|: %e dx:%e dy: %e  delta(T): %e\n\n",k7,fabs(x7),fabs(y7-cos(T)),dx7,dy7,del7);
  572.         fprintf(out,"eps 10e-9 k: %d, |x(T)|: %e  |y(T)-cos(T)|: %e dx:%e dy: %e  delta(T): %e\n\n",k9,fabs(x9),fabs(y9-cos(T)),dx9,dy9,del9);
  573.         fprintf(out,"eps 10e-11 k: %d, |x(T)|: %e  |y(T)-cos(T)|: %e dx:%e dy: %e  delta(T): %e\n\n",k11,fabs(x11),fabs(y11-cos(T)),dx11,dy11,del11);
  574.         fprintf(out,"Runge number x for T: %g\n\n", fabs((x7-x9)/(x9-x11)));
  575.         fprintf(out,"Runge number y for T: %g\n\n", fabs((y7-y9)/(y9-y11)));
  576.         R*=10;
  577.     }
  578.  
  579.     return y;
  580. }
  581.  
  582. int main(void){
  583.  
  584.     FILE*out=fopen("output.txt","w");
  585.     //pristrel(out);
  586.     zapuskrunge(out);
  587.     //test(out);
  588.    
  589.     return sqrt(9);
  590. }
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!