C語言 例題5-6 四階 Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 , 真實解 W(t)= -2e^(-t) + t ^2 - 2t + 3

C語言 例題5-6 四階 Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 , 真實解　W(t)= -2e^(-t) + t ^2 - 2t + 3


/****************************************************************/
/* Module : RK4.c                         */
/* Section: 10.2                         */
/* Cheney-Kincaid, Numerical Mathematics and Computing, 5th ed, */
/* Brooks/Cole Publ. Co.                                        */
/* Copyright (c) 2003.  All rights reserved.                    */
/* For educational use with the Cheney-Kincaid textbook.        */
/* Absolutely no warranty implied or expressed.                 */ 
/*                                 */
/* Description: Runge-Kutta 4th Order             */
/*                                 */
/****************************************************************/

#include <stdio.h>
#include <stdlib.h>

float f(float t, float y); /* define prototype */


float f(float t, float y)
{
  return (-y+t*t+1);
}


void rk4(float f(float, float), float t, float y, float h, int n)
{
  int k;
  float F1, F2, F3, F4, ta;

  ta = t;

  for (k = 1; k <= n; k++)
  {
    F1 = h * f(t,y);
    F2 = h * f(t + 0.5 * h, y + F1 * 0.5);
    F3 = h * f(t + 0.5 * h, y + F2 * 0.5);
    F4 = h * f(t + h, y + F3);
    y += (F1 + 2 * (F2 + F3) + F4) / 6.0;

    t = ta + k * h;
    if (k%10==0)
        printf("k = %d\tt = %0.2lf \ty = %0.8lf\n", k, t, y);
  }

}


void main()
{
    const int n = 100;
    const float a = 0.0, b = 1.0;
    float h,t=0, y = 1.0;

    h = (b - a) / n;
    t = a;
    printf("k = %d\tt = %0.2lf \ty = %0.8lf\n", 0, t, y);
    rk4 (f, t, y, h, n);

}


輸出畫面
k = 0 t = 0.00 y = 1.00000000
k = 10 t = 0.10 y = 1.00032520
k = 20 t = 0.20 y = 1.00253856
k = 30 t = 0.30 y = 1.00836384
k = 40 t = 0.40 y = 1.01936018
k = 50 t = 0.50 y = 1.03693891
k = 60 t = 0.60 y = 1.06237674
k = 70 t = 0.70 y = 1.09682941
k = 80 t = 0.80 y = 1.14134204
k = 90 t = 0.90 y = 1.19686067
k = 100 t = 1.00 y = 1.26424098