学习数值方法解常微分方程的笔记

工程上在满足精度要求的前提下，多用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 - y1
    dydt[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;
}