C語言 例題5-4 利用 四階 Runge-Jutta Method 求ODE y'=-y+x+1 y(0)=1 , 0<= t <=1 取 h=0.2

C語言 例題5-4 利用 四階 Runge-Jutta Method   求ODE y'=-y+x+1  y(0)=1 , 0<= t <=1    取 h=0.2

源自於

https://www.geeksforgeeks.org/runge-kutta-4th-order-method-solve-differential-equation/

Given following inputs,

• An ordinary differential equation that defines value of dy/dx in the form x and y.
• Initial value of y, i.e., y(0)

Thus we are given below.

The task is to find value of unknown function y at a given point x.

The Runge-Kutta method finds approximate value of y for a given x. Only first order ordinary differential equations can be solved by using the Runge Kutta 4th order method.

Below is the formula used to compute next value yn+1 from previous value yn. The value of n are 0, 1, 2, 3, ….(x – x0)/h. Here h is step height and xn+1 = x0 + h

 Lower step size means more accuracy.

The formula basically computes next value yn+1 using current yn plus weighted average of four increments.

• k1 is the increment based on the slope at the beginning of the interval, using y
• k2 is the increment based on the slope at the midpoint of the interval, using y + hk1/2.
• k3 is again the increment based on the slope at the midpoint, using using y + hk2/2.
• k4 is the increment based on the slope at the end of the interval, using y + hk3.

/***********************************************
****EULER METHOD FOR DIFFERENTIAL EQUATIONS*****
* C語言 例題5-4 利用 四階 Runge-Jutta Method
* 求ODE y'=-y+x+1  y(0)=1 , 0<= t <=1
* 取 h=0.2
*
***********************************************/
#include<stdio.h>
#include<math.h>

// A sample differential equation "dy/dx = (x - y)/2"
float dydx(float x, float y)
{
   return -y+x+1;
}

// Finds value of y for a given x using step size h
// and initial value y0 at x0.
float rungeKutta(float x0, float y0, float x, float h)
{
    // Count number of iterations using step size or
    // step height h
    int n = (int)((x - x0) / h);

    float k1, k2, k3, k4, k5;

    // Iterate for number of iterations
    float y = y0;
    for (int i=1; i<=n+1; i++)
    {
        // Apply Runge Kutta Formulas to find
        // next value of y
        k1 = h*dydx(x0, y);
        k2 = h*dydx(x0 + 0.5*h, y + 0.5*k1);
        k3 = h*dydx(x0 + 0.5*h, y + 0.5*k2);
        k4 = h*dydx(x0 + h, y + k3);

        // Update next value of y
        y = y + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4);;
        printf("\nThe value of y=%0.5lf at x=%0.2lf (4-order Rung-Kuttea: k1=%0.4lf ,k2=%0.4lf ,k3=%0.4lf ,k4=%0.4lf ) "
        ,y,x0, k1,k2,k3,k4);
        // Update next value of x
        x0 = x0 + h;
    }

    return y;
}

// Driver method
int main()
{
    float x0 = 0, y = 1, x =1.0, h = 0.2;
     rungeKutta(x0, y, x, h);
    return 0;
}

輸出畫面
The value of y=1.01873 at x=0.00 (4-order Rung-Kuttea: k1=0.0000 ,k2=0.0200 ,k3=0.0180 ,k4=0.0364 )
The value of y=1.07032 at x=0.20 (4-order Rung-Kuttea: k1=0.0363 ,k2=0.0526 ,k3=0.0510 ,k4=0.0661 )
The value of y=1.14882 at x=0.40 (4-order Rung-Kuttea: k1=0.0659 ,k2=0.0793 ,k3=0.0780 ,k4=0.0903 )
The value of y=1.24933 at x=0.60 (4-order Rung-Kuttea: k1=0.0902 ,k2=0.1012 ,k3=0.1001 ,k4=0.1102 )
The value of y=1.36789 at x=0.80 (4-order Rung-Kuttea: k1=0.1101 ,k2=0.1191 ,k3=0.1182 ,k4=0.1265 )
The value of y=1.50120 at x=1.00 (4-order Rung-Kuttea: k1=0.1264 ,k2=0.1338 ,k3=0.1330 ,k4=0.1398 )