C語言 例題5-6 四階 Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 求y(0.2)= ??

C語言 例題5-6 四階 Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1  求y(0.2)= ??

/*
C code using Runge-Kutta 4th order method
Problem: Here we have to find y(0.2) and y(0.4),
Given dy/dx=-y+x^2+1 where y=1 when x=0

Algorithm:
Step 1: input x0,y0,h,last point n

Step 2:m1=f(xi,yi)

Step 3:m2=f(xi+h/2,yi+m1h/2)

Step 4:m3=f(xi+h/2,yi+m2h/2)

Step 5:m4=f(xi+h,yi+m3h)

Step 6:yi+1=yi+(m1+2m2+2m3+m4/6)h

Step 5:Display output
*/

#include<stdio.h>
#include <math.h>

#define F(x,y)  -y +1 + (x*x)

int  main()
{
    double y0,x0,y1,n,h,f,k1,k2,k3,k4;
    x0=0;
    y0=1;
    h=0.01;
    n=0.2;
   
    printf("\nGiven dy/dx=-y+x^2+1 where y(0)=1 ");
    printf("\n The value of last point:%0.2lf\n ",n);

    for(; x0<n; x0=x0+h)
    {
        f=F(x0,y0);
        k1 = h * f;
        f = F(x0+h/2,y0+k1/2);
        k2 = h * f;
        f = F(x0+h/2,y0+k2/2);
        k3 = h * f;
        f = F(x0+h/2,y0+k2/2);
        k4 = h * f;
        y1 = y0 + ( k1 + 2*k2 + 2*k3 + k4)/6;
        printf("\n\n  k1 = %.4lf  ",k1);
        printf("\n\n  k2 = %.4lf ",k2);
        printf("\n\n  k3 = %.4lf ",k3);
        printf("\n\n  k4 = %.4lf ",k4);
        printf("\n\n  y(%.4lf) = %.3lf ",x0+h,y1);
        y0=y1;
    }
    return 0;
}


 輸出畫面
 Given dy/dx=-y+x^2+1 where y(0)=1
 The value of last point:0.20


  k1 = 0.0000 

  k2 = 0.0001

  k3 = 0.0001

  k4 = 0.0001

  y(0.0100) = 1.000

  k1 = 0.0000 

  k2 = 0.0002

  k3 = 0.0002

  k4 = 0.0002

  y(0.0200) = 1.000

  k1 = 0.0000 

  k2 = 0.0002

  k3 = 0.0003

  k4 = 0.0003

  y(0.0300) = 1.000

  k1 = 0.0000 

  k2 = 0.0003

  k3 = 0.0003

  k4 = 0.0003

  y(0.0400) = 1.001

  k1 = 0.0000 

  k2 = 0.0004

  k3 = 0.0004

  k4 = 0.0004

  y(0.0500) = 1.001

  k1 = 0.0000 

  k2 = 0.0005

  k3 = 0.0005

  k4 = 0.0005

  y(0.0600) = 1.001

  k1 = 0.0000 

  k2 = 0.0006

  k3 = 0.0006

  k4 = 0.0006

  y(0.0700) = 1.002

  k1 = 0.0000 

  k2 = 0.0007

  k3 = 0.0007

  k4 = 0.0007

  y(0.0800) = 1.003

  k1 = 0.0000 

  k2 = 0.0008

  k3 = 0.0008

  k4 = 0.0008

  y(0.0900) = 1.003

  k1 = 0.0000 

  k2 = 0.0009

  k3 = 0.0009

  k4 = 0.0009

  y(0.1000) = 1.004

  k1 = 0.0001 

  k2 = 0.0010

  k3 = 0.0010

  k4 = 0.0010

  y(0.1100) = 1.005

  k1 = 0.0001 

  k2 = 0.0011

  k3 = 0.0011

  k4 = 0.0011

  y(0.1200) = 1.006

  k1 = 0.0001 

  k2 = 0.0012

  k3 = 0.0012

  k4 = 0.0012

  y(0.1300) = 1.007

  k1 = 0.0001 

  k2 = 0.0013

  k3 = 0.0013

  k4 = 0.0013

  y(0.1400) = 1.008

  k1 = 0.0001 

  k2 = 0.0014

  k3 = 0.0014

  k4 = 0.0014

  y(0.1500) = 1.009

  k1 = 0.0001 

  k2 = 0.0015

  k3 = 0.0015

  k4 = 0.0015

  y(0.1600) = 1.010

  k1 = 0.0002 

  k2 = 0.0016

  k3 = 0.0016

  k4 = 0.0016

  y(0.1700) = 1.012

  k1 = 0.0002 

  k2 = 0.0016

  k3 = 0.0016

  k4 = 0.0016

  y(0.1800) = 1.013

  k1 = 0.0002 

  k2 = 0.0017

  k3 = 0.0017

  k4 = 0.0017

  y(0.1900) = 1.015

  k1 = 0.0002 

  k2 = 0.0018

  k3 = 0.0018

  k4 = 0.0018

  y(0.2000) = 1.016