C語言 例題5-11 剛性微分方程 y'(t)=-15y(t) , 0<=t , y(0)=1 真實解 y(t)= exp(-15t)

C語言   例題5-11 剛性微分方程   y'(t)=-15y(t)  , 0<=t  , y(0)=1
 真實解  y(t)= exp(-15t)


在數學領域中，剛性方程（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)=-15y(t)  , 0<=t  , y(0)=1
 // 真實解  y(t)= exp(-15t)

 #include <stdio.h>
 #include <math.h>
 #define  F(y,t)   (-15*y)
 #define  W(t)     (exp(-15*t))

 void main()
 {
    int i,n=1000;
    double a=0.0,b=1.0,y0=1.0,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.0000000  1.0000000  0.0000000
0.10  0.2231302  0.2231302  0.0000000
0.20  0.0497871  0.0497871  0.0000000
0.30  0.0111090  0.0111090  0.0000000
0.40  0.0024788  0.0024788  0.0000000
0.50  0.0005531  0.0005531  0.0000000
0.60  0.0001234  0.0001234  0.0000000
0.70  0.0000275  0.0000275  0.0000000
0.80  0.0000061  0.0000061  0.0000000
0.90  0.0000014  0.0000014  0.0000000
1.00  0.0000003  0.0000003  0.0000000

Command exited with non-zero status 233

https://www.bluffton.edu/homepages/facstaff/nesterd/java/slopefields.html