最近在网上搜了很多基于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;
}
效果显示:

在这里插入图片描述

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

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐