C語言 例題 5-7 利用 二階 Runge-Kutta 解 二階常微分程式 y'' + 4y' + 4y = 4 cos(t) + 3 sin(t) , y(0)=1 , y'(0)=0

C語言 例題 5-7 利用 二階 Runge-Kutta 解 二階常微分程式
y'' + 4y' + 4y = 4 cos(t) + 3 sin(t)   ,  y(0)=1  , y'(0)=0

設 z= y'

z' + 4z + 4y = 4 cos(t) + 3 sin(t)

z ' = f2( y , z , t) = - 4z - 4y + 4 cos(t) + 3 sin(t) 
f'  = f1(y,z ,t) = z  ,  y(0)=1

設 y1= y , y2=z

y'1 = y' = f1(y1,y2,t) = y2  ,  y1(0)=1   << f'  = f1(y,z ,t) = z  ,  y(0)=1 >>

y'2 = z' = f2(y1,y2,t) = - 4y2 - 4y1 + 4 cos(t) + 3 sin(t)  , y2(0)=0   <<- 4z - 4y + 4 cos(t) + 3 sin(t) >>


/* ex5-7.c based on Second-Order Runge-Kutta Method
 * to approximate the solution of the m order
 * system of first-order initial-value problem.
 *    y1=f1(y1,y2,...,ym,t)
 *    y2=f2(y1,y2,...,ym,t)
 *        .
 *        .
 *    ym=fm(y1,y2,...,ym,t)
 *    a<=t<=b, y1(a)=y01,y2(a)=y02,...,ym(a)=y0m.
 * at (n+1) equally spaced numbers in the interval
 * [a,b].
 */
#include <stdio.h>
#include <math.h>

#define f1(y1,y2,t)   (y2)
#define f2(y1,y2,t)   ((-4.0*y2-4.0*y1)+4*cos(t)+3*sin(t))

void main()
{
   int i=0,j,n=4,m=2;
   double t,a=0.0,b=0.4,h,k1[10],k2[10],y[10],y0[10];
   h=(b-a)/n;
   t=a;
   y0[1]=1.0;
   y0[2]=0.0;
 
   for(j=1;j<=m;j++)
      y[j]=y0[j];
     
   printf("\ni =%1d \n",i); 
   printf("t         y1         y2 \n");
   printf("=============================\n");
   printf("%.2lf %10.8lf %10.8lf \n",t,y0[1],y0[2]);
   for(i=1;i<=n;i++)
   {
       
        printf("\ni =%1d \n",i);
        for(j=1;j<=m;j++)
        {
        if(j==1)
            k1[j]=h*f1(0,y[j+1],0);
        else if(j==2)
                k1[j]=h*f2(y[j-1],y[j],t);
               
            printf("====================\n");
            printf("t         k1%1d\n",j);
            printf("%.2lf %10.8lf  \n",t,k1[j]);   
        }
     
        for(j=1;j<=m;j++)
        {
        if(j==1)
                k2[j]=h*f1(0,(y[j+1]+k1[j+1]),0);
        else if(j==2)
            k2[j]=h*f2((y[j-1]+k1[j-1]),(y[j]+k1[j]),(t+h));
           
            printf("====================\n");
            printf("t         k2%1d\n",j);
            printf("%.2lf %10.8lf  \n",t+h,k2[j]);               
           
        }

        for(j=1;j<=m;j++)
        y[j]=y[j]+0.5*(k1[j]+k2[j]);
     
        t=a+i*h;
        if(i%1==0)
            {
                printf("\nt         y1         y2 \n");
                printf("=============================\n");
                printf("%.2lf %10.8lf %10.8lf\n\n",t,y[1],y[2]);
            }
   }
   return;
}



輸出畫面
i =0
t         y1         y2
=============================
0.00 1.00000000 0.00000000

i =1
====================
t         k11
0.00 0.00000000 
====================
t         k12
0.00 0.00000000 
====================
t         k21
0.10 0.00000000 
====================
t         k22
0.10 0.02795169 

t         y1         y2
=============================
0.10 1.00000000 0.01397585


i =2
====================
t         k11
0.10 0.00139758 
====================
t         k12
0.10 0.02236135 
====================
t         k21
0.20 0.00363372 
====================
t         k22
0.20 0.03653352 

t         y1         y2
=============================
0.20 1.00251565 0.04342328


i =3
====================
t         k11
0.20 0.00434233 
====================
t         k12
0.20 0.03325186 
====================
t         k21
0.30 0.00766751 
====================
t         k22
0.30 0.03737741 

t         y1         y2
=============================
0.30 1.00852057 0.07873791


i =4
====================
t         k11
0.30 0.00787379 
====================
t         k12
0.30 0.03588726 
====================
t         k21
0.40 0.01146252 
====================
t         k22
0.40 0.03284208 

t         y1         y2
=============================
0.40 1.01818873 0.11310259