【C语言】实现 4阶(经典)龙格-库塔法 求解二阶微分方程

最近在网上搜了很多基于C语言实现龙格-库塔法 求解二阶微分方程的相关文章，发现基本都没解决我的问题，代码也没开源，于是自己找了很多相关资料写了一个具体实现的例子方便大家理解，大家***窥一斑而知全豹***

问题如下：

解题思路：

只需要将以上式子转化成以下一阶常微分方程来求即可

实现代码：

#include <stdio.h>

double function(double x , double y[] , int j){
    if(j == 1)
        return x * y[1] + y[0] ;    //如果j = 1 ,f = x * y[1] + y[0]，然后返回回去
    return y[1] ;   //否则f = y[1]，然后返回f即可
}

void rungekutta(double x ,  double *y , double h){
    double ywork[2] , k0[2] , k1[2] , k2[2] , k3[2] ;
    int j ;
    for(j = 0 ; j < 2 ; ++ j)
        k0[j] = h * function(x , y , j) ;   //计算k1
    for(j = 0 ; j < 2 ; ++ j)
        ywork[j] = y[j] + 0.5 * k0[j] ;    
         //用数组ywork存储y的变化量
    for(j = 0 ; j < 2 ; ++ j)
        k1[j] = h * function(x + 0.5 * h , ywork , j);//计算k2
    for(j = 0 ; j < 2 ; ++ j)
        ywork[j] = y[j] + 0.5 * k1[j] ;   
         //将y的变化量存储到ywork中
    for(j = 0 ; j < 2 ; ++ j)
        k2[j] = h * function(x + 0.5 * h , ywork , j);//计算k3
    for(j = 0 ; j < 2 ; ++ j)
        ywork[j] = y[j] + k2[j] ;   
        //更新ywork数组，存储y的变化量
    for(j = 0 ; j < 2 ; ++ j)
        k3[j] = h * function(x + h , ywork , j) ;   //计算k4
    for(j = 0 ; j < 2 ; ++ j)
        y[j] = y[j] + (k0[j] + 2 * k1[j] + 2 * k2[j] + k3[j]) / 6;     //计算y0和y1，用j循环先求y0，再求y1
}

int main() {
    double h , x , y[2] ; //声明变量，y0,y1使用数组y[0],y[1]表示
    h = 0.1 ;
    x = 0 ;
    y[0] = 1 , y[1] = 1 ;   //设置运算初始值
    int i; 
    for(i = 0 ; i < 10 ; ++ i){
        rungekutta(x , y , h) ; 
        //调用Runge Kutta函数进行运算，传入需要用到的数值：x , y0 , y1 , h ;
        x = x + h ;
        printf("n = %d: x = %lf, y0 = %lf, y1 = %lf\n" , i + 1 , x , y[0] , y[1]) ;
    }
    return 0;
}

效果显示：

大家需要详细流程图的可以评论我，因学习繁重这里笔者就写这么多。