c语言-四阶龙格-库塔法

#include<stdio.h>
#include<math.h>
#define n 14
//double func1(double x, double y);
double func2(double x, double y);
int main(){
    double a = 0.0, b = 1.4;//求解区间为[a,b]
    double h, m = 0.0;//h为步长，取0.1；  m为y的初值。
    double x[n + 1] = { 0.0 }, y[n + 1] = { 0.0 };//初始化数组
    h = (b - a) / n;//求步长
    x[0] = a;
    y[0] = m;
    printf("要求解的常微分方程为：dy/dx=1+y^2,初值为y(0)=0.0;
");
    printf("下面为预估校正法结果：
");
    //下面为预估校正法（即向前欧拉法和梯形法的结合，也称改进的欧拉法）
    for (int i = 1; i <= n; i++)
    {
        y[i] = y[i - 1] + h*(1 + y[i - 1] * y[i - 1]);
        x[i] = a + i*h;
        y[i] = y[i - 1] + h*0.5*((1 + y[i - 1] * y[i - 1]) + (1 + y[i] * y[i]));
        printf("x[%d]=%e , y[%d]=%e
", i, x[i], i, y[i]);
    }
    printf("下面为四阶龙格-库塔法结果：
");
    double k1 = 0.0, k2 = 0.0, k3 = 0.0, k4 = 0.0;
    for (int i = 1; i <= n; i++)
    {
        x[i] = 0.0, y[i] = 0.0;
    }
    x[0] = a;
    y[0] = m;
    for (int i = 1; i <= n; i++)
    {
        k1 = h*func2(x[i - 1], y[i - 1]);
        k2 = h*func2(x[i - 1] + h / 2.0, y[i - 1] + 0.5*k1);
        k3 = h*func2(x[i - 1] + h / 2.0, y[i - 1] + 0.5*k2);
        k4 = h*func2(x[i - 1] + h, y[i - 1] + k3);
        y[i] = y[i - 1] + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4);
        printf("x[%d]=%e , y[%d]=%e
", i, x[i], i, y[i]);
        x[i] = a + i*h;
    }
    return 0;
}
/*
double func1(double x, double y)//dy/dx=x-y+1,func1=x-y+1;
{
double f = 0;
f = x - y + 1;
return f;
}
*/
double func2(double x, double y)//dy/dx=1+y*y,func2=1+y*y;
{
    double f = 0.0;
    f = 1 + y*y;
    return f;
}

-*--------------------------

R-K法一般形式：