【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;
}
效果显示:
大家需要详细流程图的可以评论我,因学习繁重这里笔者就写这么多。
更多推荐
已为社区贡献1条内容
所有评论(0)