API tools faq
paste
Guest User

Untitled

a guest
Apr 25th, 2018
81
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 15.03 KB | None | 0 0
raw download clone embed print report
  1. #include "stdio.h"
  2. #include "conio.h"
  3. #include "math.h"
  4. #include<stdlib.h>
  5.  
  6. #define RAND_MAX    2147483647
  7. typedef double (*Function)(double x, double y);
  8. typedef double (*Function2)(double t, double x, double y);
  9. typedef double (*Function1)(double x);
  10.  
  11. double myFunc(double x, double y){
  12.     //return 1 + 0.3*y + sin(x) - y*y;
  13.     return cos(1.5*x+y)+1.5*(x-y);
  14. }
  15.  
  16.  
  17. double y1(double x, double y1, double y2){
  18.         return y2;
  19. }
  20. double y2(double x, double y1, double y2){
  21.         return (1/(x*x)-1)*y1-y2/x;
  22. }
  23.  
  24. //for me a and b == 2
  25.  
  26. //f
  27. double sysOne(double t, double x, double y){
  28.     return exp((-1)*(x*x + t*t)) + 2*x;
  29. }
  30.  
  31. //g
  32. double sysTwo(double t, double x, double y){
  33.     return 2*y*y + t;
  34. }
  35.  
  36. double funcP(double x){ return 2*x; }
  37. double funcG(double x){ return 2; }
  38. double funcF(double x){ return 4*x; }
  39.  
  40. void resultOutput(double *y,double a,double b,double h){
  41.         double x = a;
  42.         int n=(b-a)/h+0.5;
  43.         printf(" N      X[i]       Y[i]\n");
  44.         for (int i=0;i<=n;i++){
  45.                 printf("%2d %10.6f %10.6f\n",i,x,y[i]);
  46.                 x+=h;
  47.         }
  48.         printf("\n");
  49. }
  50.  
  51. void systemResultOutput(double *x, double *y, double a, double b, double h){
  52.         double t = a;
  53.         int n=(b-a)/h;
  54.         printf(" i      X[i]       Y[i]       Z[i]\n");
  55.         for (int i=0;i<=n;i++){
  56.                 printf("%2d %10.6f %10.6f %10.6f\n",i,t,x[i],y[i]);
  57.                 t+=h;
  58.         }
  59.         printf("\n");
  60. }
  61.  
  62. double rungeMethod(Function f, double a, double b, double y0, double h, double eps, double *y)
  63. {
  64.     double delta;
  65.     int n = (b-a)/h, i;
  66.     y[0] = y0;
  67.     double *tmp = new double[n+1];//temp array
  68.     int multiply = 1; //mul
  69.     do {
  70.     //points values
  71.     double y1 = y0;
  72.     for (i=1;i<=n;i++){
  73.         for (int j = 0; j < multiply; j++){
  74.         double x = a + ((i-1)*multiply+j)*h;
  75.         //koeff
  76.         double k1 = h * f(x,y1);
  77.         double k2 = h * f(x+h/2,y1+k1/2);
  78.         double k3 = h * f(x+h/2,y1+k2/2);
  79.         double k4 = h * f(x+h,y1+k3);
  80.         double y2 = y1 + (k1+2*k2+2*k3+k4)/6;
  81.         y1 = y2;
  82.         }
  83.         tmp[i] = y1;
  84.     }
  85.     //checks eps
  86.     if (multiply == 1){
  87.         delta = eps+1;
  88.     }else{
  89.         delta = 0;
  90.         for (i=1;i<=n;i++){
  91.         delta += fabs(tmp[i]-y[i]);
  92.         }
  93.     }
  94.     //transfer to new array
  95.     if (delta>eps){
  96.         for (i=1; i<=n; i++) y[i]=tmp[i];
  97.         multiply *= 2; //mul
  98.         h/=2;
  99.     }
  100.     } while (delta>eps);
  101.     return h;
  102. }
  103.  
  104. double methodRKSystem(Function2 f, Function2 g, double a, double b,
  105.                     double x0, double y0, double h, double eps, double *x, double *y)
  106. {
  107.     double delta;
  108.     int n = (b-a)/h, i;
  109.     x[0] = x0;
  110.     y[0] = y0;
  111.     double *tmpx = new double[n+1]; //tmp
  112.     double *tmpy = new double[n+1]; //tmp
  113.     int multiply = 1;
  114.     do {
  115.     //point values
  116.     double x1 = x0;
  117.     double y1 = y0;
  118.     for (i=1;i<=n;i++){
  119.         for (int j = 0; j < multiply; j++){
  120.         double t = a + ((i-1)*multiply+j)*h;
  121.         //koeff
  122.         double k1 = h * f(t,x1,y1);
  123.         double l1 = h * g(t,x1,y1);
  124.         double k2 = h * f(t+h/2,x1+k1/2,y1+l1/2);
  125.         double l2 = h * g(t+h/2,x1+k1/2,y1+l1/2);
  126.         double k3 = h * f(t+h/2,x1+k2/2,y1+l2/2);
  127.         double l3 = h * g(t+h/2,x1+k2/2,y1+l2/2);
  128.         double k4 = h * f(t+h,x1+k3,y1+l3);
  129.         double l4 = h * g(t+h,x1+k3,y1+l3);
  130.         //new values
  131.         double x2 = x1 + (k1+2*k2+2*k3+k4)/6;
  132.         double y2 = y1 + (l1+2*l2+2*l3+l4)/6;
  133.         x1 = x2;
  134.         y1 = y2;
  135.         }
  136.         tmpx[i] = x1;
  137.         tmpy[i] = y1;
  138.     }
  139.     //eps
  140.     if (multiply==1){
  141.  
  142.         delta = eps+1;
  143.     }else{
  144.         delta = 0;
  145.         for (i=1;i<=n;i++){
  146.         delta += fabs(tmpx[i]-x[i]);
  147.         delta += fabs(tmpy[i]-y[i]);
  148.         }
  149.     }
  150.     //new array
  151.     if (delta>eps){
  152.         for (i=1;i<=n;i++) { x[i]=tmpx[i]; y[i]=tmpy[i]; }
  153.         multiply *= 2;
  154.         h/=2;
  155.     }
  156.     } while (delta>eps);
  157.     return h;
  158. }
  159.  
  160.  
  161. double milnMethod(Function f, double a, double b, double y0, double h, double eps, double *y){
  162.         double delta, tmp[5], dy[5];
  163.         int div = 1;
  164.         int n = (b-a)/h, i, j, k;
  165.         y[0] = y0;
  166.         double *prev = new double[n+1];
  167.         do {
  168.                 //4 points by runge
  169.                 rungeMethod(f, a, a+h*3, y0, h, eps, tmp);
  170.                 double x = a;
  171.                 for (i=0; i<=n; i++) {
  172.                         y[i] = tmp[0];
  173.                         for (k=0; k<div; k++){
  174.                                 //derivatives
  175.                                 for(j=0; j<4; j++) dy[j] = f(x+h*j,tmp[j]);
  176.                                 //extrapolation
  177.                                 double yn = tmp[0] + 4*h/3*(2*dy[3]-dy[2]+2*dy[1]);
  178.                                 //interpolation
  179.                                 dy[4] = f(x+h*4,yn);
  180.                                 yn = tmp[2] + h/3*(dy[2]+4*dy[3]+dy[4]);
  181.                                 //1 step right
  182.                                 tmp[4] = yn;
  183.                                 x += h;
  184.                                 for (int i=0;i<4;i++) tmp[i] = tmp[i+1];
  185.                         }
  186.                 }
  187.                 //check
  188.                 if (div==1){
  189.                         delta = eps+1;
  190.                 }else{
  191.                         delta = 0;
  192.                         for (i=0;i<n;i++) delta += fabs(y[i]-prev[i]);
  193.                 }
  194.                 if (delta>eps){
  195.                         div *= 2;
  196.                         h /= 2;
  197.                         for (i=0;i<n;i++) prev[i]=y[i];
  198.                 }
  199.         } while (delta>eps);
  200.         return h;
  201. }
  202.  
  203. double adamsMethod(Function f, double a, double b, double y0, double h, double eps, double *y){
  204.         double delta, tmp[5], dy[5];
  205.         int div = 1;
  206.         int n = (b-a)/h, i, j, k;
  207.         y[0] = y0;
  208.         double *prev = new double[n+1];
  209.         do {
  210.                 //4 points frm Runge
  211.                 rungeMethod(f, a, a+h*3, y0, h, eps, tmp);
  212.                 double x = a;
  213.                 for (i=0; i<=n; i++) {
  214.                         y[i] = tmp[0];
  215.                         for (k=0; k<div; k++){
  216.                                 //derivatives
  217.                                 for(j=0; j<4; j++) dy[j] = f(x+h*j, tmp[j]);
  218.                                 //extrapolation
  219.                                 double yn = tmp[3] + h/24*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0]);
  220.                                 //interpolation
  221.                                 dy[4] = f(x+h*4,yn);
  222.                                 yn = tmp[3] + h/24*(9*dy[4]+19*dy[3]-5*dy[2]+dy[1]);
  223.                                 //1 step right
  224.                                 tmp[4] = yn;
  225.                                 x += h;
  226.                                 for(int i=0; i<4; i++) tmp[i] = tmp[i+1];
  227.                         }
  228.                 }
  229.                 //check
  230.                 if (div==1){
  231.                         delta = eps+1;
  232.                 }else{
  233.                         delta = 0;
  234.                         for (i=0;i<n;i++) delta += fabs(y[i]-prev[i]);
  235.                 }
  236.                 if (delta>eps){
  237.                         //h to the next
  238.                         div*=2;
  239.                         h/=2;
  240.                         for (i=0;i<n;i++) prev[i]=y[i];
  241.                 }
  242.         } while (delta>eps);
  243.         return h;
  244. }
  245.  
  246.  
  247. void runMethod(Function1 p, Function1 g, Function1 f, double a, double b, double h,
  248.                 double a0, double a1, double a2, double b0, double b1, double b2, double *y)
  249. {
  250.         int n = ((b-a)/h+0.5)+1, i, j ;
  251.         static double c[100],d[100];
  252.         //straight
  253.         double m0 = -2+h*p(a);
  254.         double n0 = 1-h*p(a)+h*h*g(a);
  255.         c[0] = (a1-a0*h)/(m0*(a1-a0*h)+n0*a1);
  256.         d[0] = n0*a2*h/(a1-a0*h)+f(a)*h*h;
  257.         for (i=1;i<n;i++){
  258.                 m0 = -2+h*p(a+h*i);
  259.                 n0 = 1-h*p(a+h*i)+h*h*g(a+h*i);
  260.                 c[i] = 1/(m0-n0*c[i-1]);
  261.                 d[i] = f(a+h*i)*h*h-n0*c[i-1]*d[i-1];
  262.         }
  263.         //feedback
  264.         y[n-1] = (b1*c[n-3]*d[n-3]+b2*h)/(b1*(1+c[n-3])+b0*h);
  265.         for (i=n-2;i>0;i--){
  266.                 y[i] = c[i-1]*(d[i-1]-y[i+1]);
  267.         }
  268.         y[0] = (a1*y[1]-a2*h)/(a1-a0*h);
  269. }
  270.  
  271. void MKR(Function1 p, Function1 g, Function1 f, double a, double b, double h,
  272.                 double a0, double a1, double a2, double b0, double b1, double b2, double *y)
  273. {
  274.         //matrix & right part vector
  275.         int n = ((b-a)/h+0.5)+1, i, j ;
  276.         static double m[50][50];
  277.         static double v[100];
  278.         for (i=0;i<n;i++){
  279.                 for (j=0;j<n;j++) m[i][j]=0;
  280.                 v[i]=0;
  281.         }
  282.         //main YR
  283.         for (i=1;i<n-1;i++){
  284.                 m[i][i-1] = (1/h-p(a+i*h)/2)/h;
  285.                 m[i][i] = -2/(h*h)+g(a+i*h);
  286.                 m[i][i+1] = 1/(h*h)+p(a+i*h)/(2*h);
  287.                 v[i] = f(a+i*h);
  288.         }
  289.         //edge conditions
  290.         m[0][0] = a0-a1/h; m[0][1] = a1/h; v[0] = a2;
  291.         m[n-1][n-1] = b0+b1/h; m[n-1][n-2] = -b1/h; v[n-1] = b2;
  292.         //solving system
  293.         for (i=0;i<n;i++){
  294.                 //div str at diag el
  295.                 double diag = m[i][i];
  296.                 for (j=0;j<n;j++) m[i][j]/=diag;
  297.                 m[i][i]=1; v[i]/=diag;
  298.                 //rem i column by substr
  299.                 for (j=i+1;j<n;j++){
  300.                         double kf = m[j][i];
  301.                         for (int k=i;k<n;k++) m[j][k]-=kf*m[i][k];
  302.                         m[j][i]=0;
  303.                         v[j]-=kf*v[i];
  304.                 }
  305.         }
  306.         //feedback
  307.         for (i=n-1;i>=0;i--){
  308.                 double yy = v[i];
  309.                 for (j=i+1;j<n;j++) yy-=y[j]*m[i][j];
  310.                 y[i] = yy;
  311.         }
  312. }
  313.  
  314.  
  315. double Compare(double *y, double *y2, double a, double b, double h){
  316.         int n = (b-a)/h+0.5;
  317.         double delta = 0;
  318.         for (int i=0;i<n;i++){
  319.                 double d = fabs(y[i]-y2[i*2]);
  320.                 if (d>delta) delta = d;
  321.         }
  322.         return delta;
  323. }
  324.  
  325. int shootMethod(Function2 y1, Function2 y2, double a, double b, double h, double eps, double g1,
  326.           double g2, double g3, double g4, double g5, double g6, double *y1r, double *y2r){
  327.         int n = (b-a)/h+0.5;
  328.         randomize();
  329.         //double L1 = (qrand()*1.0)/RAND_MAX, L2 = (rand()*1.0)/RAND_MAX;
  330.         double L1 = (rand()*1.0)/RAND_MAX, L2 =10+ (rand()*1.0)/RAND_MAX;
  331.         double psi1, psi2;
  332.         int i=1;
  333.         //val at 1st point
  334.         if (g2!=0){
  335.                 methodRKSystem(y1, y2, a, b, L1, (g3-g1*L1)/g2, h, eps, y1r, y2r);
  336.         }else{
  337.                 methodRKSystem(y1, y2, a, b, (g3-g2*L1)/g1, L1, h, eps, y1r, y2r);
  338.         }
  339.         psi1 = g4*y1r[n-1]+g5*y2r[n-1]-g6;
  340.         while (1){
  341.                 //koshi y1(a)=L
  342.                 if (g2!=0){
  343.                         methodRKSystem(y1, y2, a, b, L2,(g3-g1*L2)/g2, h, eps, y1r, y2r);
  344.                 }else{
  345.                         methodRKSystem(y1, y2, a, b, (g3-g2*L2)/g1, L2, h, eps, y1r, y2r);
  346.                 }
  347.                 psi2 = g4*y1r[n-1]+g5*y2r[n-1]-g6;
  348.                 if (fabs(psi1)+fabs(psi2)<eps) break;
  349.                 //corrected L
  350.                 double L = L2 - (L2-L1)*psi2/(psi2-psi1);
  351.                 L1 = L2; L2 = L;
  352.                 psi1 = psi2;
  353.                 i++;
  354.         }
  355.         return i;
  356. }
  357.  
  358. void output(double *x, double *y, double a, double b, double h){
  359.         double t = a;
  360.         int n=(b-a)/h;
  361.         printf(" i      X[i]       Y1[i]      Y2[i]\n");
  362.         for (int i=0;i<=n;i++){
  363.                 printf("%2d %10.6f %10.6f %10.6f\n",i,t,x[i],y[i]);
  364.                 t+=h;
  365.         }
  366.         printf("\n");
  367. }
  368.  
  369. int main()
  370. {
  371.     int result = 0,key=0;
  372.     double a, b, h, eps, y0, x0;
  373.     double a0, a1, a2, b0, b1, b2;
  374.     double y[100], x[100];
  375.     double g1, g2, g3, g4, g5, g6;
  376.  
  377.     while(1)
  378.     {
  379.        clrscr();
  380.        printf("\n1) - ЊҐв®¤ ђг­ЈҐ-Љгвв .");
  381.        printf("\n2) - ЊҐв®¤ ђг­ЈҐ-Љгвв  ¤«п бЁб⥬ га ў­Ґ­Ё©.");
  382.        printf("\n3) - ЊҐв®¤ Ђ¤ ¬б .");
  383.        printf("\n4) - ЊҐв®¤ ЊЁ«­ .");
  384.        printf("\n5) - ЊҐв®¤ Є®­Ґз­ле а §­®б⥩.");
  385.        printf("\n6) - ЊҐв®¤ Їа®Ј®­ЄЁ.");
  386.        printf("\n7) - ЊҐв®¤ бв५мЎл.");
  387.        printf("\n0) - ‚л室");
  388.        printf("\nЌ ¦¬ЁвҐ Є« ўЁ¦г ўлЎа ­­®Ј® ў ¬Ё Їг­Єв  ¬Ґ­о. ");
  389.        //int menuOpt = 0;
  390.        //scanf("%d", &menuOpt);
  391.        key = getch();
  392.        switch (key)
  393.        {
  394.             case 48: exit(1);
  395.             case 49:
  396.             {
  397.                printf("\n\n‚ўҐ¤ЁвҐ a, b, eps, h, y(0): ");
  398.                scanf("%lg %lg %lg %lg %lg", &a, &b, &eps, &h, &y0);
  399.                rungeMethod(myFunc, a, b, y0, h, eps, y);
  400.                resultOutput(y, a, b, h);
  401.                break;
  402.             }
  403.             case 50:
  404.             {
  405.                printf("\n\n‚ўҐ¤ЁвҐ a, b, eps, h, y(0), x(0): ");
  406.                scanf("%lg %lg %lg %lg %lg %lg", &a, &b, &eps, &h, &y0, &x0);
  407.                methodRKSystem(sysOne, sysTwo, a, b, x0, y0, h, eps, x, y);
  408.                systemResultOutput(x, y, a, b, h);
  409.                break;
  410.             }
  411.            case 51:
  412.            {
  413.               printf("\n\n‚ўҐ¤ЁвҐ a, b, eps, h, y(0): ");
  414.               scanf("%lg %lg %lg %lg %lg", &a, &b, &eps, &h, &y0);
  415.               adamsMethod(myFunc, a, b, y0, h, eps, y);
  416.               resultOutput(y, a, b, h);
  417.               break;
  418.            }
  419.            case 52:
  420.            {
  421.  
  422.                printf("\n\n‚ўҐ¤ЁвҐ a, b, eps, h, y(0): ");
  423.                scanf("%lg %lg %lg %lg %lg", &a, &b, &eps, &h, &y0);
  424.                milnMethod(myFunc, a, b, y0, h, eps, y);
  425.                resultOutput(y, a, b, h);
  426.                break;
  427.            }
  428.            case 53:
  429.            {
  430.               printf("\n\n‚ўҐ¤ЁвҐ a, b, h: ");
  431.               scanf("%lg %lg %lg", &a, &b, &h);
  432.               eps=0.000001; a0=1; a1=-0.5; a2=2; b0=0; b1=1; b2=4;
  433.               MKR(funcP, funcG, funcF, a, b, h, a0, a1, a2, b0, b1, b2, y);
  434.               MKR(funcP, funcG, funcF, a, b, h/2, a0, a1, a2, b0, b1, b2, x);
  435.               resultOutput(y, a, b, h);
  436.               printf("\nЏ®ЈаҐи­®бвм: %10.6f\n", Compare(y,x,a,b,h)/3);
  437.               break;
  438.            }
  439.            case 54:
  440.            {
  441.               printf("\n\n‚ўҐ¤ЁвҐ a, b: ");
  442.               scanf("%lg %lg", &a, &b);
  443.               eps=0.001; h=0.05; a0=1; a1=-0.5; a2=2; b0=0; b1=1; b2=4;
  444.               runMethod(funcP, funcG, funcF, a, b, h, a0, a1, a2, b0, b1, b2, y);
  445.               runMethod(funcP, funcG, funcF, a, b, h/2, a0, a1, a2, b0, b1, b2, x);
  446.               resultOutput(y, a, b, h);
  447.               printf("\nЏ®ЈаҐи­®бвм: %10.6f\n", Compare(y,x,a,b,h)/3);
  448.               break;
  449.            }
  450.            case 55:
  451.            {
  452.                printf("\n\n‚ўҐ¤ЁвҐ a, b, h: ");
  453.                scanf("%lg %lg %lg", &a, &b, &h);
  454.                eps=0.000001; g1=0; g2=1; g3=1.2; g4=1; g5=0; g6=2.1;
  455.                shootMethod(y1, y2, a, b, h, eps, g1, g2, g3, g4, g5, g6, x, y);
  456.                output(x, y, a, b, h);
  457.                break;
  458.            }
  459.         }
  460.        printf("\n\nЌ ¦¬ЁвҐ «оЎго Є« ўЁиг...");
  461.        getch();
  462.     }
  463.    return 0;
  464. }
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!