本文主要是介绍学习数值方法解常微分方程的笔记，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
工程上在满足精度要求的前提下，多用4阶龙格库塔法解常微分方程
取步长0.001，步数10，C语言代码如下
#include <stdio.h>
#include <math.h>// 定义微分方程系统的导数函数 f(t, y)
double f(double t, double y) {return y;
}// 四阶龙格-库塔方法（RK4）
void rk4(double (*f)(double, double), double t0, double y0, double h, int n_steps, double *results) {int i;double k1, k2, k3, k4;results[0] = y0; // 存储初始值for (i = 1; i <= n_steps; ++i) {// 计算四个阶段的斜率k1 = f(t0 + (i - 1) * h, results[i - 1]);k2 = f(t0 + (i - 1) * h + 0.5 * h, results[i - 1] + 0.5 * h * k1);k3 = f(t0 + (i - 1) * h + 0.5 * h, results[i - 1] + 0.5 * h * k2);k4 = f(t0 + (i - 1) * h + h, results[i - 1] + h * k3);// 使用四阶龙格-库塔公式更新结果results[i] = results[i - 1] + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);}
}int main() {const double t0 = 0.0; // 初始时间const double y0 = 1.0; // 初始状态const double h = 0.001; // 时间步长const int n_steps = 10; // 总迭代步数double results[n_steps + 1]; // 用于存储结果的时间序列rk4(f, t0, y0, h, n_steps, results);for(int i=0;i<n_steps;i++){printf("%.9lf %.9lf\n",results[i],results[i]-exp(i*h));}return 0;
}结果如下
1.000000000 0.000000000
1.001000500 0.000000000
1.002002001 0.000000000
1.003004505 0.000000000
1.004008011 0.000000000
1.005012521 0.000000000
1.006018036 0.000000000
1.007024557 0.000000000
1.008032086 0.000000000
1.009040622 0.000000000与精确值的误差几乎没有
对于多阶常微分方程，可以逐级化为一阶的
比如y''' + 2y'' + 2y' + y = 80，化为y1' = y2, y2' = y3, y3' = 80 - 2y3 - 2y2 - y1，这个80 - 2y3 - 2y2 - y1相当于一阶里的f(t,y)
#include <stdio.h>
#include <math.h>// 定义三阶常微分方程的导数函数 f(t, y, y', y'', y''')
void f(double t, double *y, double *dydt) {// 假设我们有一个三阶常微分方程：y''' + 2y'' + 2y' + y = 80// 将其转化为一阶系统：y1' = y2, y2' = y3, y3' = 80 - 2y3 - 2y2 - y1dydt[0] = y[1];dydt[1] = y[2];dydt[2] = 80 - 2 * y[2] - 2 * y[1] - y[0];
}// 四阶龙格-库塔方法（RK4）针对一阶系统
void rk4_system(void (*f)(double, double*, double*), double t0, double *y0, double h, int n_steps, double *results) {int i, j;double k1[3], k2[3], k3[3], k4[3];for (j = 0; j < 3; ++j) {results[j] = y0[j]; // 存储初始值}for (i = 1; i <= n_steps; ++i) {// 计算四个阶段的斜率f(t0 + (i - 1) * h, &results[i * 3 - 3], k1);// 计算中间状态（临时变量）double y_mid[3];for (j = 0; j < 3; ++j) {y_mid[j] = results[i * 3 - 3 + j] + 0.5 * h * k1[j]; // 修正：逐元素乘法并计算中间状态}f(t0 + (i - 1) * h + 0.5 * h, y_mid, k2);f(t0 + (i - 1) * h + 0.5 * h, y_mid, k3); // 修正：使用中间状态计算斜率，此处重复调用f函数，需要删除其中一个// 计算下一个时间步的中间状态double next_y_mid[3];for (j = 0; j < 3; ++j) {next_y_mid[j] = results[i * 3 - 3 + j] + h * k3[j]; // 修正：逐元素乘法并计算下一个时间步的中间状态}f(t0 + (i - 1) * h + h, next_y_mid, k4); // 修正：使用下一个时间步的中间状态计算斜率// 使用四阶龙格-库塔公式更新结果for (j = 0; j < 3; ++j) {results[i * 3 + j] = results[i * 3 - 3 + j] + (h / 6.0) * (k1[j] + 2.0 * k2[j] + 2.0 * k3[j] + k4[j]);}}
}int main() {const double t0 = 0.0; // 初始时间double y0[] = {0.0, 0.0, 0.0}; // 初始状态const double h = 0.1; // 时间步长const int n_steps = 100; // 总迭代步数double results[3 * (n_steps + 1)]; // 用于存储结果的时间序列rk4_system(f, t0, y0, h, n_steps, results);for(int i=0;i<n_steps;i+=10){printf("%lf",results[3*i]);printf("\n");}return 0;
}这篇关于学习数值方法解常微分方程的笔记的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
