C語言 例題5-11 y'(t)=(50.0/y)-50*y , 0<=t , y(0)= sqrt(2) 真實解 y(t)= sqrt(1+(1.0/exp(100*t)))

 C語言 例題5-11   y'(t)=(50.0/y)-50*y    , 0<=t  , y(0)= sqrt(2)
                              真實解  y(t)=  sqrt(1+(1.0/exp(100*t)))

/* ex5-11.c based on Four-Order Runge-Kutta
 * Method to approximate the solution of the
 * initial-value problem
 *   y'=f(y,t), a<=t<=b, y(a)=y0
 * at (n+1) equally spaced numbers in the interval
 * [a,b]: input a,b,n,and initial condition y0.
 */
 //  y'(t)=(50.0/y)-50*y    , 0<=t  , y(0)= sqrt(2)
 // 真實解  y(t)=  sqrt(1+(1.0/exp(100*t)))

 #include <stdio.h>
 #include <math.h>
 #define  F(y,t)   ((50.0/y)-50*y)
 #define  W(t)     (sqrt(1+(1.0/exp(100*t))))

 void main()
 {
    int i,n=10;
    double a=0.0,b=1.0,y0=sqrt(2),k1,k2,k3,k4,h,t,y,err;
    h=(b-a)/n;
    t=a;
    y=y0;
    err=fabs(y-W(t));
    printf("t      y(t)       w(t)      error\n");
    printf("=====================================\n");
    printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
    for(i=1;i<=n;i++)
    {
       k1=h*F(y,t);
       k2=h*F((y+k1/2.0),(t+h/2.0));
       k3=h*F((y+k2/2.0),(t+h/2.0));
       k4=h*F((y+k3),(t+h));
       y=y+(k1+2*k2+2*k3+k4)/6.0;
     
       t=a+i*h;
       err=fabs(y-W(t));
       if(i%1==0)
          printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
    }
    return;
}


輸出畫面
t               y(t)                              w(t)                                            error
==================================================
0.00  1.4142136                        1.4142136                             0.0000000
0.10 -15.8525966                      1.0000227                           16.8526193
0.20 -215.7458639                    1.0000000                         216.7458639
0.30 -2957.4012711                  1.0000000                        2958.4012711
0.40 -40541.0340389                1.000000                        40542.0340389
0.50 -555750.0076712              1.0000000                     555751.0076712
0.60 -7618406.3551152            1.0000000                   7618407.3551152
0.70 -104435653.7847015        1.0000000               104435654.7847015
0.80 -1431638753.9652824      1.0000000             1431638754.9652824
0.90 -19625381252.2740898    1.0000000           19625381253.2740898
1.00 -269031267999.9240723  1.0000000         269031268000.9240723

Command exited with non-zero status 11


剛性方程（stiffness equation）是指一個微分方程，其數值分析的解只有在時間間隔很小時才會穩定，只要時間間隔略大，其解就會不穩定。


/* ex5-11.c based on Four-Order Runge-Kutta
 * Method to approximate the solution of the
 * initial-value problem
 *   y'=f(y,t), a<=t<=b, y(a)=y0
 * at (n+1) equally spaced numbers in the interval
 * [a,b]: input a,b,n,and initial condition y0.
 */
 //  y'(t)=(50.0/y)-50*y    , 0<=t  , y(0)= sqrt(2)
 // 真實解  y(t)=  sqrt(1+(1.0/exp(100*t)))

 #include <stdio.h>
 #include <math.h>
 #define  F(y,t)   ((50.0/y)-50*y)
 #define  W(t)     (sqrt(1+(1.0/exp(100*t))))

 void main()
 {
    int i,n=1000;
    double a=0.0,b=1.0,y0=sqrt(2),k1,k2,k3,k4,h,t,y,err;
    h=(b-a)/n;
    t=a;
    y=y0;
    err=fabs(y-W(t));
    printf("t      y(t)       w(t)      error\n");
    printf("=====================================\n");
    printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
    for(i=1;i<=n;i++)
    {
       k1=h*F(y,t);
       k2=h*F((y+k1/2.0),(t+h/2.0));
       k3=h*F((y+k2/2.0),(t+h/2.0));
       k4=h*F((y+k3),(t+h));
       y=y+(k1+2*k2+2*k3+k4)/6.0;
     
       t=a+i*h;
       err=fabs(y-W(t));
       if(i%100==0)
          printf("%.2lf %10.7lf %10.7lf %10.7lf\n" ,t,y,W(t),err);
    }
    return;
}

輸出畫面
t      y(t)       w(t)      error
=====================================
0.00  1.4142136  1.4142136  0.0000000
0.10  1.0000227  1.0000227  0.0000000
0.20  1.0000000  1.0000000  0.0000000
0.30  1.0000000  1.0000000  0.0000000
0.40  1.0000000  1.0000000  0.0000000
0.50  1.0000000  1.0000000  0.0000000
0.60  1.0000000  1.0000000  0.0000000
0.70  1.0000000  1.0000000  0.0000000
0.80  1.0000000  1.0000000  0.0000000
0.90  1.0000000  1.0000000  0.0000000
1.00  1.0000000  1.0000000  0.0000000

Command exited with non-zero status 233

     h=0.1

h=0.01