嗯.

记得那是2021年4月6日(周二, 清明节假期), 15:22. 我在吃饭. 有个师妹问我怎么用一个叫 GSL 的 C++ library 来解常微分方程组(ODEs). 我查了一下,

The GNU Scientific Library (GSL) is a numerical library for C and C++ programmers. It is free software under the GNU General Public License.

https://www.gnu.org/software/gsl/

按 GSL 说明文档, GSL 是一个跟 scipy 类似的, 可用作科学计算的库. 任何科学计算的库的基本功能都是类似的, 所以这对我来说不是一个很大的问题. 奈何我有拖延症, 遇到师妹更是如此. 于是直到三个星期过后, 我才开始细看说明文档, 动笔写这篇文章.

闲话休题.

前置工作

我的操作系统是 windows 7, 前置工作大体有以下四个:

  • 把 C++ 安装在 windows 系统;
  • 把 GSL 安装在 windows 系统;
  • 安装 visual studio code, 配置 C++ 的编程环境;
  • 配置 visual studio code, 使之能正常调用 GSL .

这个过程说复杂其实也不复杂. 有机会出个 windows 版的安装教程.

这过程折腾了我6个小时... 从晚12点到晨6点. 其中编译安装 GSL 大概花了两个小时. 期间师妹多次喊我去睡觉, 我拒绝了.

待解 ODEs

跟 MATLAB 或 mathematica 相比, GSL 的 ODEs 求解过程有自己的风格. GSL 要求 ODEs 要写成如下的形式:

$$\frac{\text{d}y_i}{\text{d}t} = f_i(t, y_1,\cdots,y_n), $$

其中 $i=1,\cdots,n$. 初值条件为

$$y_1(t_0) = a_1, \cdots, y_n(t_0) = a_n.$$

求解过程中某些算法会用到 Jacobi 矩阵

$$J_{ij} = \frac{\partial f_i}{\partial y_j}(t, y_1,\cdots,y_n),$$

但并不是所有算法都需要. 故下文不对 Jacobi 矩阵的相关内容展开说明.

为了行文方便, 设待解的 ODEs 是

$$\frac{\text{d}y_1}{\text{d}t} = -y_2, \frac{\text{d}y_2}{\text{d}t} = y_1.$$

初值条件是

$$y_1(0) = 1, y_1(0) = 0.$$

显然, 该 ODEs 的唯一解是

$$y_1 = \cos t, y_2 = \sin t.$$

限制 $t\in[0, 2\pi]$. 最后可以检验数值解的准确性.

求解 ODEs (复杂版)

定义函数 $f_i$ 和矩阵 $J_{ij}$

GSL 对函数 $f_i$ 和矩阵 $J_{ij}$ 的封装形式有特定的要求.

函数 $f_i$ 的定义:

// 函数自变量是 y, 因变量储存在 f, 函数返回 GSL_SUCCESS (=0) 表示计算成功.
int func (double t, const double y[], double f[], void *params) {
    (void)(t); // 假装用过变量 t. 缺少这部可能会在编译时报 warning.
    double mu = *(double *)params; // 可通过 params 导入函数所需参数, 非必要. 
    f[0] = - mu*y[1];
    f[1] =   mu*y[0];
    return GSL_SUCCESS; // GSL规定, 函数运行成功返回 GSL_SUCCESS, 类似返回 0.
}

Jacobi 矩阵 $J_{ij}$ 的定义:

// 某些步进函数需要用到 Jacobi 矩阵, 本文不需要用到.
int jac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
    double mu = *(double *)params;
    gsl_matrix *m = gsl_matrix_alloc (2, 2);
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 1, -1.0*mu);
    gsl_matrix_set (m, 1, 0, 1.0*mu);
    gsl_matrix_set (m, 1, 1, 0.0);
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;
    return GSL_SUCCESS;
}

定义 ODEs 系统

ODEs 系统 (ODEs system) 是个由函数 $f_i$ , Jacobi 矩阵 $J_{ij}$, 方程维数 (dim) 及相关参数 (params) 构成的四元组 (4-tuple). 其中 Jacobi 矩阵和相关参数不是必需的.

// 定义 ODEs 系统
double mu = 1;
size_t dim = 2;
// gsl_odeiv2_system sys = {func, jac, 2, &mu}; // 本文不需要 Jacobi 矩阵
gsl_odeiv2_system sys = {func, NULL, dim, &mu};

定义步进函数

步进函数 (Stepping function) 确定解 ODEs 所涉及的算法. GSL 中所涉及的算法基本都是 Runge-Kutta 法, 不同的算法精度可能差异很大. 步进函数的相关说明如下:

步进函数是最底层的部件, 它能使得解从 $t$ 步进到 $t+h$, 其中 $h$ (相对)固定, 并估计结果的局部误差.

The lowest level components are the stepping functions which advance a solution from time $t$ to $t+h$ for a fixed step-size h and estimate the resulting local error.

步长 $h$ 在计算过程中给定, 这里不需要给定. 相关代码:

const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk2;
gsl_odeiv2_step * step = gsl_odeiv2_step_alloc (T, dim);
// gsl_odeiv2_step_rk8pd 表示求解 ODEs 所用的算法, 可以换成 gsl_odeiv2_step_type 类中别的算法, 其中有些算法要求 Jacobi 矩阵.

定义控制函数

控制函数 (control function) 是用来调控解的误差的, 使得解的误差落在一定范围内. 不同的控制函数控制不同类型的误差. 其的相关说明如下:

步进函数生成 (ODEs 的) 解, 而 (步长) 控制函数检验这些解(是否符合)给定的范围, 并尝试按用户给定的误差来确定最优的步长.

The control function examines the proposed change to the solution produced by a stepping function and attempts to determine the optimal step-size for a user-specified level of error.

相关代码:

// 定义控制函数, 总误差 = 绝对误差 + 相对误差
double eps_abs = 1e-3; // 绝对误差
double eps_rel = 0; // 相对误差
gsl_odeiv2_control * con = gsl_odeiv2_control_y_new (eps_abs, eps_rel);

定义演进函数

演进函数 (evolution function 进化函数? 这个该怎么翻译...) 是用来控制解的步进. 相关说明如下:

演化函数通过结合步进函数和控制函数, 以及给定的步长, 让解步进一步.

The evolution function combines the results of a stepping function and control function to reliably advance the solution forward one step using an acceptable step-size.

相关代码如下:

gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (dim);
// dim 表示 ODEs 的维度.
//按说明文档, gsl_odeiv2_evolve 类似乎只有一个对象 gsl_odeiv2_evolve_alloc .

这步是计算前的最后准备, 相当于告诉计算机演进所涉及的 ODEs 是 $2$ 维的.

笔者不明白这一步的意义何在... 仅仅是告诉计算机方程的维数吗?

步进计算

有了上述准备, 便可以开始计算了.

由于 C++ 无法直接返回数组, 因此 GSL 不像其他的科学计算库那样, 给定方程和条件就可以直接返回结果. GSL 只能实现用已知 ODEs 在 $t$ 处的值(初值), 来计算 $t+h$ 处的值. 显然 $h$ 的值不能过大或者过小. 要求 ODEs 在特定区间上的数值解, 还需要在外层套个循环.

相关代码

int status = gsl_odeiv2_evolve_apply_fixed_step(e, con, step, &sys, &t, h, y);
// status == 0 表示计算成功
// 计算结果储存在变量 t 和 y

完整代码和返回结果

完整代码如下:

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <math.h>

int func (double t, const double y[], double f[], void *params) {
    (void)(t); // 没有用到 t, 避免某些 warning.
    double mu = *(double *)params;
    f[0] = - mu*y[1];
    f[1] =   mu*y[0];
    return GSL_SUCCESS;
}

void print_one_step(int n, double t, double y[]){
    // 打印计算结果
    printf("step %d: t, y1, y2\n", n);
    printf ("pred: %1.5f, %1.5f, %1.5f\n",   t, y[0],   y[1]);   // 预测值
    printf ("real: %1.5f, %1.5f, %1.5f\n\n", t, cos(t), sin(t)); // 真实值
}

int main(void){

    // 定义 ODEs 系统
    double mu = 1;
    size_t dim = 2;
    gsl_odeiv2_system sys = {func, NULL, dim, &mu};

    //定义步进函数
    const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk2;
    gsl_odeiv2_step * step = gsl_odeiv2_step_alloc (T, dim);

    // 定义控制函数, 总误差 = 绝对误差 + 相对误差
    double eps_abs = 1e-3; // 绝对误差
    double eps_rel = 0; // 相对误差
    gsl_odeiv2_control * con = gsl_odeiv2_control_y_new (eps_abs, eps_rel);

    // 定义演进函数
    gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (dim);

    // 开始计算
    double t  = 0.0 ; // 初值
    double t1 = 2*M_PI; // 终止值
    double h  = M_PI / 40; // 步长
    double y[2] = { 1.0, 0.0 }; // 初值

    int n = 0;
    print_one_step(n, t, y);

    while (t < t1) {
        int status = gsl_odeiv2_evolve_apply_fixed_step(e, con, step, &sys, &t, h, y);
        if (status != GSL_SUCCESS){
            // 异常则退出计算
            printf("error, status = %d\n", status);
            break;
        }
        n ++;
        print_one_step(n, t, y);
    }

    // 释放内存
    gsl_odeiv2_step_free(step);
    gsl_odeiv2_control_free(con);
    gsl_odeiv2_evolve_free(e);

    return 0;
}

返回结果如下:

step 0: t, y1, y2
pred: 0.00000, 1.00000, 0.00000
real: 0.00000, 1.00000, 0.00000

step 1: t, y1, y2
pred: 0.07854, 0.99692, 0.07846
real: 0.07854, 0.99692, 0.07846

step 2: t, y1, y2
pred: 0.15708, 0.98769, 0.15643
real: 0.15708, 0.98769, 0.15643

step 3: t, y1, y2
pred: 0.23562, 0.97237, 0.23344
real: 0.23562, 0.97237, 0.23345

step 4: t, y1, y2
pred: 0.31416, 0.95105, 0.30902
real: 0.31416, 0.95106, 0.30902

step 5: t, y1, y2
pred: 0.39270, 0.92387, 0.38268
real: 0.39270, 0.92388, 0.38268

step 6: t, y1, y2
pred: 0.47124, 0.89100, 0.45399
real: 0.47124, 0.89101, 0.45399

step 7: t, y1, y2
pred: 0.54978, 0.85263, 0.52249
real: 0.54978, 0.85264, 0.52250

step 8: t, y1, y2
pred: 0.62832, 0.80901, 0.58778
real: 0.62832, 0.80902, 0.58779

step 9: t, y1, y2
pred: 0.70686, 0.76039, 0.64944
real: 0.70686, 0.76041, 0.64945

step 10: t, y1, y2
pred: 0.78540, 0.70709, 0.70710
real: 0.78540, 0.70711, 0.70711

step 11: t, y1, y2
pred: 0.86394, 0.64944, 0.76039
real: 0.86394, 0.64945, 0.76041

step 12: t, y1, y2
pred: 0.94248, 0.58777, 0.80900
real: 0.94248, 0.58779, 0.80902

step 13: t, y1, y2
pred: 1.02102, 0.52249, 0.85262
real: 1.02102, 0.52250, 0.85264

step 14: t, y1, y2
pred: 1.09956, 0.45398, 0.89099
real: 1.09956, 0.45399, 0.89101

step 15: t, y1, y2
pred: 1.17810, 0.38267, 0.92386
real: 1.17810, 0.38268, 0.92388

step 16: t, y1, y2
pred: 1.25664, 0.30901, 0.95103
real: 1.25664, 0.30902, 0.95106

step 17: t, y1, y2
pred: 1.33518, 0.23344, 0.97234
real: 1.33518, 0.23345, 0.97237

step 18: t, y1, y2
pred: 1.41372, 0.15643, 0.98766
real: 1.41372, 0.15643, 0.98769

step 19: t, y1, y2
pred: 1.49226, 0.07845, 0.99689
real: 1.49226, 0.07846, 0.99692

step 20: t, y1, y2
pred: 1.57080, -0.00000, 0.99997
real: 1.57080, 0.00000, 1.00000

step 21: t, y1, y2
pred: 1.64934, -0.07846, 0.99688
real: 1.64934, -0.07846, 0.99692

step 22: t, y1, y2
pred: 1.72788, -0.15643, 0.98765
real: 1.72788, -0.15643, 0.98769

step 23: t, y1, y2
pred: 1.80642, -0.23344, 0.97233
real: 1.80642, -0.23345, 0.97237

step 24: t, y1, y2
pred: 1.88496, -0.30901, 0.95102
real: 1.88496, -0.30902, 0.95106

step 25: t, y1, y2
pred: 1.96350, -0.38267, 0.92384
real: 1.96350, -0.38268, 0.92388

step 26: t, y1, y2
pred: 2.04204, -0.45397, 0.89097
real: 2.04204, -0.45399, 0.89101

step 27: t, y1, y2
pred: 2.12058, -0.52248, 0.85260
real: 2.12058, -0.52250, 0.85264

step 28: t, y1, y2
pred: 2.19911, -0.58776, 0.80898
real: 2.19911, -0.58779, 0.80902

step 29: t, y1, y2
pred: 2.27765, -0.64942, 0.76037
real: 2.27765, -0.64945, 0.76041

step 30: t, y1, y2
pred: 2.35619, -0.70708, 0.70707
real: 2.35619, -0.70711, 0.70711

step 31: t, y1, y2
pred: 2.43473, -0.76037, 0.64941
real: 2.43473, -0.76041, 0.64945

step 32: t, y1, y2
pred: 2.51327, -0.80898, 0.58775
real: 2.51327, -0.80902, 0.58779

step 33: t, y1, y2
pred: 2.59181, -0.85260, 0.52247
real: 2.59181, -0.85264, 0.52250

step 34: t, y1, y2
pred: 2.67035, -0.89096, 0.45396
real: 2.67035, -0.89101, 0.45399

step 35: t, y1, y2
pred: 2.74889, -0.92383, 0.38266
real: 2.74889, -0.92388, 0.38268

step 36: t, y1, y2
pred: 2.82743, -0.95100, 0.30900
real: 2.82743, -0.95106, 0.30902

step 37: t, y1, y2
pred: 2.90597, -0.97231, 0.23343
real: 2.90597, -0.97237, 0.23345

step 38: t, y1, y2
pred: 2.98451, -0.98763, 0.15642
real: 2.98451, -0.98769, 0.15643

step 39: t, y1, y2
pred: 3.06305, -0.99686, 0.07845
real: 3.06305, -0.99692, 0.07846

step 40: t, y1, y2
pred: 3.14159, -0.99994, -0.00000
real: 3.14159, -1.00000, 0.00000

step 41: t, y1, y2
pred: 3.22013, -0.99685, -0.07846
real: 3.22013, -0.99692, -0.07846

step 42: t, y1, y2
pred: 3.29867, -0.98762, -0.15643
real: 3.29867, -0.98769, -0.15643

step 43: t, y1, y2
pred: 3.37721, -0.97230, -0.23343
real: 3.37721, -0.97237, -0.23345

step 44: t, y1, y2
pred: 3.45575, -0.95099, -0.30900
real: 3.45575, -0.95106, -0.30902

step 45: t, y1, y2
pred: 3.53429, -0.92381, -0.38266
real: 3.53429, -0.92388, -0.38268

step 46: t, y1, y2
pred: 3.61283, -0.89094, -0.45396
real: 3.61283, -0.89101, -0.45399

step 47: t, y1, y2
pred: 3.69137, -0.85257, -0.52246
real: 3.69137, -0.85264, -0.52250

step 48: t, y1, y2
pred: 3.76991, -0.80895, -0.58774
real: 3.76991, -0.80902, -0.58779

step 49: t, y1, y2
pred: 3.84845, -0.76034, -0.64940
real: 3.84845, -0.76041, -0.64945

step 50: t, y1, y2
pred: 3.92699, -0.70705, -0.70705
real: 3.92699, -0.70711, -0.70711

step 51: t, y1, y2
pred: 4.00553, -0.64939, -0.76035
real: 4.00553, -0.64945, -0.76041

step 52: t, y1, y2
pred: 4.08407, -0.58773, -0.80895
real: 4.08407, -0.58779, -0.80902

step 53: t, y1, y2
pred: 4.16261, -0.52245, -0.85257
real: 4.16261, -0.52250, -0.85264

step 54: t, y1, y2
pred: 4.24115, -0.45395, -0.89093
real: 4.24115, -0.45399, -0.89101

step 55: t, y1, y2
pred: 4.31969, -0.38265, -0.92380
real: 4.31969, -0.38268, -0.92388

step 56: t, y1, y2
pred: 4.39823, -0.30898, -0.95097
real: 4.39823, -0.30902, -0.95106

step 57: t, y1, y2
pred: 4.47677, -0.23342, -0.97228
real: 4.47677, -0.23345, -0.97237

step 58: t, y1, y2
pred: 4.55531, -0.15641, -0.98760
real: 4.55531, -0.15643, -0.98769

step 59: t, y1, y2
pred: 4.63385, -0.07845, -0.99682
real: 4.63385, -0.07846, -0.99692

step 60: t, y1, y2
pred: 4.71239, 0.00001, -0.99991
real: 4.71239, 0.00000, -1.00000

step 61: t, y1, y2
pred: 4.79093, 0.07846, -0.99682
real: 4.79093, 0.07846, -0.99692

step 62: t, y1, y2
pred: 4.86947, 0.15643, -0.98759
real: 4.86947, 0.15643, -0.98769

step 63: t, y1, y2
pred: 4.94801, 0.23343, -0.97227
real: 4.94801, 0.23345, -0.97237

step 64: t, y1, y2
pred: 5.02655, 0.30899, -0.95096
real: 5.02655, 0.30902, -0.95106

step 65: t, y1, y2
pred: 5.10509, 0.38265, -0.92378
real: 5.10509, 0.38268, -0.92388

step 66: t, y1, y2
pred: 5.18363, 0.45395, -0.89091
real: 5.18363, 0.45399, -0.89101

step 67: t, y1, y2
pred: 5.26217, 0.52245, -0.85255
real: 5.26217, 0.52250, -0.85264

step 68: t, y1, y2
pred: 5.34071, 0.58773, -0.80893
real: 5.34071, 0.58779, -0.80902

step 69: t, y1, y2
pred: 5.41925, 0.64938, -0.76032
real: 5.41925, 0.64945, -0.76041

step 70: t, y1, y2
pred: 5.49779, 0.70703, -0.70702
real: 5.49779, 0.70711, -0.70711

step 71: t, y1, y2
pred: 5.57633, 0.76033, -0.64937
real: 5.57633, 0.76041, -0.64945

step 72: t, y1, y2
pred: 5.65487, 0.80893, -0.58771
real: 5.65487, 0.80902, -0.58779

step 73: t, y1, y2
pred: 5.73341, 0.85255, -0.52243
real: 5.73341, 0.85264, -0.52250

step 74: t, y1, y2
pred: 5.81195, 0.89091, -0.45393
real: 5.81195, 0.89101, -0.45399

step 75: t, y1, y2
pred: 5.89049, 0.92377, -0.38263
real: 5.89049, 0.92388, -0.38268

step 76: t, y1, y2
pred: 5.96903, 0.95094, -0.30897
real: 5.96903, 0.95106, -0.30902

step 77: t, y1, y2
pred: 6.04757, 0.97225, -0.23341
real: 6.04757, 0.97237, -0.23345

step 78: t, y1, y2
pred: 6.12611, 0.98757, -0.15641
real: 6.12611, 0.98769, -0.15643

step 79: t, y1, y2
pred: 6.20465, 0.99679, -0.07844
real: 6.20465, 0.99692, -0.07846

step 80: t, y1, y2
pred: 6.28319, 0.99987, 0.00001
real: 6.28319, 1.00000, 0.00000

求解 ODEs (简化版)

驱动器 (driver) 可以把上述的步进函数, 控制函数, 演进函数打包, 用起来更加方便. 下面利用驱动器的写法, 把上述例子重写一遍.

定义驱动器

相关代码如下:

const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk2; // 算法
double h  = M_PI / 40; // 步长
double eps_abs = 1e-3; // 绝对误差
double eps_rel = 0; // 相对误差
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); // 定义驱动器, sys 是 ODEs 系统.

完整代码和返回结果

完整代码如下:

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <math.h>

int func (double t, const double y[], double f[], void *params) {
    (void)(t); 
    double mu = *(double *)params;
    f[0] = - mu*y[1];
    f[1] =   mu*y[0];
    return GSL_SUCCESS;
}

void print_one_step(int n, double t, double y[]){
    // 打印计算结果
    printf("step %d: t, y1, y2\n", n);
    printf ("pred: %1.5f, %1.5f, %1.5f\n",   t, y[0],   y[1]);   // 预测值
    printf ("real: %1.5f, %1.5f, %1.5f\n\n", t, cos(t), sin(t)); // 真实值
}

int main(void){

    // 定义 ODEs 系统
    double mu = 1;
    size_t dim = 2;
    gsl_odeiv2_system sys = {func, NULL, dim, &mu};

    const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk2; // 算法
    double h  = M_PI / 40; // 步长
    double eps_abs = 1e-3; // 绝对误差
    double eps_rel = 0; // 相对误差
    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);

    double t = 0;
    double t1 = 2*M_PI;
    double y[] = {1, 0};

    int n = 0;
    print_one_step(n, t, y);

    while (t < t1) {
        int status = gsl_odeiv2_driver_apply(d, &t, t + h, y);
        if (status != GSL_SUCCESS){
            // 异常则退出计算
            printf("error, status = %d\n", status);
            break;
        }

        n++;
        print_one_step(n, t, y);
    }

    gsl_odeiv2_driver_free(d);

    return 0;
}

返回结果如下:

step 0: t, y1, y2
pred: 0.00000, 1.00000, 0.00000
real: 0.00000, 1.00000, 0.00000

step 1: t, y1, y2
pred: 0.07854, 0.99692, 0.07846
real: 0.07854, 0.99692, 0.07846

step 2: t, y1, y2
pred: 0.15708, 0.98769, 0.15643
real: 0.15708, 0.98769, 0.15643

step 3: t, y1, y2
pred: 0.23562, 0.97237, 0.23344
real: 0.23562, 0.97237, 0.23345

step 4: t, y1, y2
pred: 0.31416, 0.95105, 0.30902
real: 0.31416, 0.95106, 0.30902

step 5: t, y1, y2
pred: 0.39270, 0.92387, 0.38268
real: 0.39270, 0.92388, 0.38268

step 6: t, y1, y2
pred: 0.47124, 0.89100, 0.45399
real: 0.47124, 0.89101, 0.45399

step 7: t, y1, y2
pred: 0.54978, 0.85263, 0.52249
real: 0.54978, 0.85264, 0.52250

step 8: t, y1, y2
pred: 0.62832, 0.80901, 0.58778
real: 0.62832, 0.80902, 0.58779

step 9: t, y1, y2
pred: 0.70686, 0.76039, 0.64944
real: 0.70686, 0.76041, 0.64945

step 10: t, y1, y2
pred: 0.78540, 0.70709, 0.70710
real: 0.78540, 0.70711, 0.70711

step 11: t, y1, y2
pred: 0.86394, 0.64944, 0.76039
real: 0.86394, 0.64945, 0.76041

step 12: t, y1, y2
pred: 0.94248, 0.58777, 0.80900
real: 0.94248, 0.58779, 0.80902

step 13: t, y1, y2
pred: 1.02102, 0.52249, 0.85262
real: 1.02102, 0.52250, 0.85264

step 14: t, y1, y2
pred: 1.09956, 0.45398, 0.89099
real: 1.09956, 0.45399, 0.89101

step 15: t, y1, y2
pred: 1.17810, 0.38267, 0.92386
real: 1.17810, 0.38268, 0.92388

step 16: t, y1, y2
pred: 1.25664, 0.30901, 0.95103
real: 1.25664, 0.30902, 0.95106

step 17: t, y1, y2
pred: 1.33518, 0.23344, 0.97234
real: 1.33518, 0.23345, 0.97237

step 18: t, y1, y2
pred: 1.41372, 0.15643, 0.98766
real: 1.41372, 0.15643, 0.98769

step 19: t, y1, y2
pred: 1.49226, 0.07845, 0.99689
real: 1.49226, 0.07846, 0.99692

step 20: t, y1, y2
pred: 1.57080, -0.00000, 0.99997
real: 1.57080, 0.00000, 1.00000

step 21: t, y1, y2
pred: 1.64934, -0.07846, 0.99688
real: 1.64934, -0.07846, 0.99692

step 22: t, y1, y2
pred: 1.72788, -0.15643, 0.98765
real: 1.72788, -0.15643, 0.98769

step 23: t, y1, y2
pred: 1.80642, -0.23344, 0.97233
real: 1.80642, -0.23345, 0.97237

step 24: t, y1, y2
pred: 1.88496, -0.30901, 0.95102
real: 1.88496, -0.30902, 0.95106

step 25: t, y1, y2
pred: 1.96350, -0.38267, 0.92384
real: 1.96350, -0.38268, 0.92388

step 26: t, y1, y2
pred: 2.04204, -0.45397, 0.89097
real: 2.04204, -0.45399, 0.89101

step 27: t, y1, y2
pred: 2.12058, -0.52248, 0.85260
real: 2.12058, -0.52250, 0.85264

step 28: t, y1, y2
pred: 2.19911, -0.58776, 0.80898
real: 2.19911, -0.58779, 0.80902

step 29: t, y1, y2
pred: 2.27765, -0.64942, 0.76037
real: 2.27765, -0.64945, 0.76041

step 30: t, y1, y2
pred: 2.35619, -0.70708, 0.70707
real: 2.35619, -0.70711, 0.70711

step 31: t, y1, y2
pred: 2.43473, -0.76037, 0.64941
real: 2.43473, -0.76041, 0.64945

step 32: t, y1, y2
pred: 2.51327, -0.80898, 0.58775
real: 2.51327, -0.80902, 0.58779

step 33: t, y1, y2
pred: 2.59181, -0.85260, 0.52247
real: 2.59181, -0.85264, 0.52250

step 34: t, y1, y2
pred: 2.67035, -0.89096, 0.45396
real: 2.67035, -0.89101, 0.45399

step 35: t, y1, y2
pred: 2.74889, -0.92383, 0.38266
real: 2.74889, -0.92388, 0.38268

step 36: t, y1, y2
pred: 2.82743, -0.95100, 0.30900
real: 2.82743, -0.95106, 0.30902

step 37: t, y1, y2
pred: 2.90597, -0.97231, 0.23343
real: 2.90597, -0.97237, 0.23345

step 38: t, y1, y2
pred: 2.98451, -0.98763, 0.15642
real: 2.98451, -0.98769, 0.15643

step 39: t, y1, y2
pred: 3.06305, -0.99686, 0.07845
real: 3.06305, -0.99692, 0.07846

step 40: t, y1, y2
pred: 3.14159, -0.99994, -0.00000
real: 3.14159, -1.00000, 0.00000

step 41: t, y1, y2
pred: 3.22013, -0.99685, -0.07846
real: 3.22013, -0.99692, -0.07846

step 42: t, y1, y2
pred: 3.29867, -0.98762, -0.15643
real: 3.29867, -0.98769, -0.15643

step 43: t, y1, y2
pred: 3.37721, -0.97230, -0.23343
real: 3.37721, -0.97237, -0.23345

step 44: t, y1, y2
pred: 3.45575, -0.95099, -0.30900
real: 3.45575, -0.95106, -0.30902

step 45: t, y1, y2
pred: 3.53429, -0.92381, -0.38266
real: 3.53429, -0.92388, -0.38268

step 46: t, y1, y2
pred: 3.61283, -0.89094, -0.45396
real: 3.61283, -0.89101, -0.45399

step 47: t, y1, y2
pred: 3.69137, -0.85257, -0.52246
real: 3.69137, -0.85264, -0.52250

step 48: t, y1, y2
pred: 3.76991, -0.80895, -0.58774
real: 3.76991, -0.80902, -0.58779

step 49: t, y1, y2
pred: 3.84845, -0.76034, -0.64940
real: 3.84845, -0.76041, -0.64945

step 50: t, y1, y2
pred: 3.92699, -0.70705, -0.70705
real: 3.92699, -0.70711, -0.70711

step 51: t, y1, y2
pred: 4.00553, -0.64939, -0.76035
real: 4.00553, -0.64945, -0.76041

step 52: t, y1, y2
pred: 4.08407, -0.58773, -0.80895
real: 4.08407, -0.58779, -0.80902

step 53: t, y1, y2
pred: 4.16261, -0.52245, -0.85257
real: 4.16261, -0.52250, -0.85264

step 54: t, y1, y2
pred: 4.24115, -0.45395, -0.89093
real: 4.24115, -0.45399, -0.89101

step 55: t, y1, y2
pred: 4.31969, -0.38265, -0.92380
real: 4.31969, -0.38268, -0.92388

step 56: t, y1, y2
pred: 4.39823, -0.30898, -0.95097
real: 4.39823, -0.30902, -0.95106

step 57: t, y1, y2
pred: 4.47677, -0.23342, -0.97228
real: 4.47677, -0.23345, -0.97237

step 58: t, y1, y2
pred: 4.55531, -0.15641, -0.98760
real: 4.55531, -0.15643, -0.98769

step 59: t, y1, y2
pred: 4.63385, -0.07845, -0.99682
real: 4.63385, -0.07846, -0.99692

step 60: t, y1, y2
pred: 4.71239, 0.00001, -0.99991
real: 4.71239, 0.00000, -1.00000

step 61: t, y1, y2
pred: 4.79093, 0.07846, -0.99682
real: 4.79093, 0.07846, -0.99692

step 62: t, y1, y2
pred: 4.86947, 0.15643, -0.98759
real: 4.86947, 0.15643, -0.98769

step 63: t, y1, y2
pred: 4.94801, 0.23343, -0.97227
real: 4.94801, 0.23345, -0.97237

step 64: t, y1, y2
pred: 5.02655, 0.30899, -0.95096
real: 5.02655, 0.30902, -0.95106

step 65: t, y1, y2
pred: 5.10509, 0.38265, -0.92378
real: 5.10509, 0.38268, -0.92388

step 66: t, y1, y2
pred: 5.18363, 0.45395, -0.89091
real: 5.18363, 0.45399, -0.89101

step 67: t, y1, y2
pred: 5.26217, 0.52245, -0.85255
real: 5.26217, 0.52250, -0.85264

step 68: t, y1, y2
pred: 5.34071, 0.58773, -0.80893
real: 5.34071, 0.58779, -0.80902

step 69: t, y1, y2
pred: 5.41925, 0.64938, -0.76032
real: 5.41925, 0.64945, -0.76041

step 70: t, y1, y2
pred: 5.49779, 0.70703, -0.70702
real: 5.49779, 0.70711, -0.70711

step 71: t, y1, y2
pred: 5.57633, 0.76033, -0.64937
real: 5.57633, 0.76041, -0.64945

step 72: t, y1, y2
pred: 5.65487, 0.80893, -0.58771
real: 5.65487, 0.80902, -0.58779

step 73: t, y1, y2
pred: 5.73341, 0.85255, -0.52243
real: 5.73341, 0.85264, -0.52250

step 74: t, y1, y2
pred: 5.81195, 0.89091, -0.45393
real: 5.81195, 0.89101, -0.45399

step 75: t, y1, y2
pred: 5.89049, 0.92377, -0.38263
real: 5.89049, 0.92388, -0.38268

step 76: t, y1, y2
pred: 5.96903, 0.95094, -0.30897
real: 5.96903, 0.95106, -0.30902

step 77: t, y1, y2
pred: 6.04757, 0.97225, -0.23341
real: 6.04757, 0.97237, -0.23345

step 78: t, y1, y2
pred: 6.12611, 0.98757, -0.15641
real: 6.12611, 0.98769, -0.15643

step 79: t, y1, y2
pred: 6.20465, 0.99679, -0.07844
real: 6.20465, 0.99692, -0.07846

step 80: t, y1, y2
pred: 6.28319, 0.99987, 0.00001
real: 6.28319, 1.00000, 0.00000

提高解的精度

笔者发现, 只要把上述的 gsl_odeiv2_step_rk2 算法换成 gsl_odeiv2_step_rk8pd, 其余不变, 就能达到很高的精度. 笔者同样不清楚其中的原理.

完整代码

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <math.h>

int func (double t, const double y[], double f[], void *params) {
    (void)(t);
    double mu = *(double *)params;
    f[0] = - mu*y[1];
    f[1] =   mu*y[0];
    return GSL_SUCCESS;
}

void print_one_step(int n, double t, double y[]){
    // 打印计算结果
    printf("step %d: t, y1, y2\n", n);
    printf ("pred: %1.19f, %1.19f, %1.19f\n",   t, y[0],   y[1]);   // 预测值
    printf ("real: %1.19f, %1.19f, %1.19f\n\n", t, cos(t), sin(t)); // 真实值
}

int main(void){

    // 定义 ODEs 系统
    double mu = 1;
    size_t dim = 2;
    gsl_odeiv2_system sys = {func, NULL, dim, &mu};

    const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk8pd; // 算法
    double h  = M_PI / 40; // 步长
    double eps_abs = 1e-3; // 绝对误差
    double eps_rel = 0; // 相对误差
    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);

    double t = 0;
    double t1 = 2*M_PI;
    double y[] = {1, 0};

    int n = 0;
    print_one_step(n, t, y);

    while (t < t1) {
        int status = gsl_odeiv2_driver_apply(d, &t, t + h, y);
        if (status != GSL_SUCCESS){
            // 异常则退出计算
            printf("error, status = %d\n", status);
            break;
        }

        n++;
        print_one_step(n, t, y);
    }

    gsl_odeiv2_driver_free(d);

    return 0;
}

返回结果

step 0: t, y1, y2
pred: 0.0000000000000000000, 1.0000000000000000000, 0.0000000000000000000
real: 0.0000000000000000000, 1.0000000000000000000, 0.0000000000000000000

step 1: t, y1, y2
pred: 0.0785398163397448280, 0.9969173337331279600, 0.0784590957278449440
real: 0.0785398163397448280, 0.9969173337331279600, 0.0784590957278449440

step 2: t, y1, y2
pred: 0.1570796326794896600, 0.9876883405951376600, 0.1564344650402308400
real: 0.1570796326794896600, 0.9876883405951377700, 0.1564344650402308700

step 3: t, y1, y2
pred: 0.2356194490192344800, 0.9723699203976765600, 0.2334453638559053600
real: 0.2356194490192344800, 0.9723699203976765600, 0.2334453638559053900

step 4: t, y1, y2
pred: 0.3141592653589793100, 0.9510565162951535300, 0.3090169943749473400
real: 0.3141592653589793100, 0.9510565162951535300, 0.3090169943749474000

step 5: t, y1, y2
pred: 0.3926990816987241400, 0.9238795325112867400, 0.3826834323650896700
real: 0.3926990816987241400, 0.9238795325112867400, 0.3826834323650897800

step 6: t, y1, y2
pred: 0.4712388980384689700, 0.8910065241883679000, 0.4539904997395466400
real: 0.4712388980384689700, 0.8910065241883679000, 0.4539904997395467500

step 7: t, y1, y2
pred: 0.5497787143782138000, 0.8526401643540922900, 0.5224985647159486900
real: 0.5497787143782138000, 0.8526401643540921800, 0.5224985647159488000

step 8: t, y1, y2
pred: 0.6283185307179586200, 0.8090169943749474500, 0.5877852522924729200
real: 0.6283185307179586200, 0.8090169943749474500, 0.5877852522924731400

step 9: t, y1, y2
pred: 0.7068583470577034500, 0.7604059656000309300, 0.6494480483301834400
real: 0.7068583470577034500, 0.7604059656000309300, 0.6494480483301836600

step 10: t, y1, y2
pred: 0.7853981633974482800, 0.7071067811865475700, 0.7071067811865473500
real: 0.7853981633974482800, 0.7071067811865475700, 0.7071067811865474600

step 11: t, y1, y2
pred: 0.8639379797371931100, 0.6494480483301837700, 0.7604059656000307100
real: 0.8639379797371931100, 0.6494480483301836600, 0.7604059656000309300

step 12: t, y1, y2
pred: 0.9424777960769379300, 0.5877852522924732500, 0.8090169943749472300
real: 0.9424777960769379300, 0.5877852522924731400, 0.8090169943749474500

step 13: t, y1, y2
pred: 1.0210176124166828000, 0.5224985647159490200, 0.8526401643540920700
real: 1.0210176124166828000, 0.5224985647159489100, 0.8526401643540921800

step 14: t, y1, y2
pred: 1.0995574287564276000, 0.4539904997395469700, 0.8910065241883676800
real: 1.0995574287564276000, 0.4539904997395468000, 0.8910065241883677900

step 15: t, y1, y2
pred: 1.1780972450961724000, 0.3826834323650899500, 0.9238795325112866300
real: 1.1780972450961724000, 0.3826834323650898400, 0.9238795325112867400

step 16: t, y1, y2
pred: 1.2566370614359172000, 0.3090169943749476200, 0.9510565162951534200
real: 1.2566370614359172000, 0.3090169943749474500, 0.9510565162951535300

step 17: t, y1, y2
pred: 1.3351768777756621000, 0.2334453638559056400, 0.9723699203976764500
real: 1.3351768777756621000, 0.2334453638559054700, 0.9723699203976765600

step 18: t, y1, y2
pred: 1.4137166941154069000, 0.1564344650402311200, 0.9876883405951375500
real: 1.4137166941154069000, 0.1564344650402309200, 0.9876883405951377700

step 19: t, y1, y2
pred: 1.4922565104551517000, 0.0784590957278452210, 0.9969173337331278500
real: 1.4922565104551517000, 0.0784590957278449990, 0.9969173337331279600

step 20: t, y1, y2
pred: 1.5707963267948966000, 0.0000000000000002914, 0.9999999999999998900
real: 1.5707963267948966000, 0.0000000000000000612, 1.0000000000000000000

step 21: t, y1, y2
pred: 1.6493361431346414000, -0.0784590957278446380, 0.9969173337331278500
real: 1.6493361431346414000, -0.0784590957278448740, 0.9969173337331279600

step 22: t, y1, y2
pred: 1.7278759594743862000, -0.1564344650402305400, 0.9876883405951376600
real: 1.7278759594743862000, -0.1564344650402308100, 0.9876883405951377700

step 23: t, y1, y2
pred: 1.8064157758141310000, -0.2334453638559050900, 0.9723699203976765600
real: 1.8064157758141310000, -0.2334453638559053400, 0.9723699203976766700

step 24: t, y1, y2
pred: 1.8849555921538759000, -0.3090169943749471200, 0.9510565162951535300
real: 1.8849555921538759000, -0.3090169943749473400, 0.9510565162951536400

step 25: t, y1, y2
pred: 1.9634954084936207000, -0.3826834323650894500, 0.9238795325112867400
real: 1.9634954084936207000, -0.3826834323650897300, 0.9238795325112867400

step 26: t, y1, y2
pred: 2.0420352248333655000, -0.4539904997395464700, 0.8910065241883679000
real: 2.0420352248333655000, -0.4539904997395466900, 0.8910065241883679000

step 27: t, y1, y2
pred: 2.1205750411731104000, -0.5224985647159485800, 0.8526401643540922900
real: 2.1205750411731104000, -0.5224985647159488000, 0.8526401643540922900

step 28: t, y1, y2
pred: 2.1991148575128552000, -0.5877852522924728000, 0.8090169943749475600
real: 2.1991148575128552000, -0.5877852522924730300, 0.8090169943749474500

step 29: t, y1, y2
pred: 2.2776546738526000000, -0.6494480483301833300, 0.7604059656000311500
real: 2.2776546738526000000, -0.6494480483301835500, 0.7604059656000310400

step 30: t, y1, y2
pred: 2.3561944901923448000, -0.7071067811865472400, 0.7071067811865477900
real: 2.3561944901923448000, -0.7071067811865474600, 0.7071067811865475700

step 31: t, y1, y2
pred: 2.4347343065320897000, -0.7604059656000306000, 0.6494480483301839900
real: 2.4347343065320897000, -0.7604059656000309300, 0.6494480483301837700

step 32: t, y1, y2
pred: 2.5132741228718345000, -0.8090169943749471200, 0.5877852522924734700
real: 2.5132741228718345000, -0.8090169943749473400, 0.5877852522924732500

step 33: t, y1, y2
pred: 2.5918139392115793000, -0.8526401643540919600, 0.5224985647159492400
real: 2.5918139392115793000, -0.8526401643540921800, 0.5224985647159489100

step 34: t, y1, y2
pred: 2.6703537555513241000, -0.8910065241883675700, 0.4539904997395471900
real: 2.6703537555513241000, -0.8910065241883677900, 0.4539904997395468600

step 35: t, y1, y2
pred: 2.7488935718910690000, -0.9238795325112865200, 0.3826834323650902300
real: 2.7488935718910690000, -0.9238795325112867400, 0.3826834323650898900

step 36: t, y1, y2
pred: 2.8274333882308138000, -0.9510565162951534200, 0.3090169943749479000
real: 2.8274333882308138000, -0.9510565162951535300, 0.3090169943749475100

step 37: t, y1, y2
pred: 2.9059732045705586000, -0.9723699203976764500, 0.2334453638559058900
real: 2.9059732045705586000, -0.9723699203976765600, 0.2334453638559055300

step 38: t, y1, y2
pred: 2.9845130209103035000, -0.9876883405951375500, 0.1564344650402313700
real: 2.9845130209103035000, -0.9876883405951376600, 0.1564344650402309800

step 39: t, y1, y2
pred: 3.0630528372500483000, -0.9969173337331278500, 0.0784590957278454710
real: 3.0630528372500483000, -0.9969173337331279600, 0.0784590957278450680

step 40: t, y1, y2
pred: 3.1415926535897931000, -0.9999999999999998900, 0.0000000000000005412
real: 3.1415926535897931000, -1.0000000000000000000, 0.0000000000000001225

step 41: t, y1, y2
pred: 3.2201324699295379000, -0.9969173337331278500, -0.0784590957278443880
real: 3.2201324699295379000, -0.9969173337331279600, -0.0784590957278448190

step 42: t, y1, y2
pred: 3.2986722862692828000, -0.9876883405951376600, -0.1564344650402302900
real: 3.2986722862692828000, -0.9876883405951377700, -0.1564344650402307300

step 43: t, y1, y2
pred: 3.3772121026090276000, -0.9723699203976765600, -0.2334453638559048100
real: 3.3772121026090276000, -0.9723699203976766700, -0.2334453638559052800

step 44: t, y1, y2
pred: 3.4557519189487724000, -0.9510565162951535300, -0.3090169943749467900
real: 3.4557519189487724000, -0.9510565162951536400, -0.3090169943749472800

step 45: t, y1, y2
pred: 3.5342917352885173000, -0.9238795325112867400, -0.3826834323650891200
real: 3.5342917352885173000, -0.9238795325112868500, -0.3826834323650896700

step 46: t, y1, y2
pred: 3.6128315516282621000, -0.8910065241883679000, -0.4539904997395461400
real: 3.6128315516282621000, -0.8910065241883679000, -0.4539904997395466900

step 47: t, y1, y2
pred: 3.6913713679680069000, -0.8526401643540922900, -0.5224985647159482400
real: 3.6913713679680069000, -0.8526401643540922900, -0.5224985647159486900

step 48: t, y1, y2
pred: 3.7699111843077517000, -0.8090169943749475600, -0.5877852522924724700
real: 3.7699111843077517000, -0.8090169943749475600, -0.5877852522924730300

step 49: t, y1, y2
pred: 3.8484510006474966000, -0.7604059656000311500, -0.6494480483301829900
real: 3.8484510006474966000, -0.7604059656000310400, -0.6494480483301835500

step 50: t, y1, y2
pred: 3.9269908169872414000, -0.7071067811865477900, -0.7071067811865469100
real: 3.9269908169872414000, -0.7071067811865476800, -0.7071067811865474600

step 51: t, y1, y2
pred: 4.0055306333269858000, -0.6494480483301843200, -0.7604059656000300400
real: 4.0055306333269858000, -0.6494480483301841000, -0.7604059656000306000

step 52: t, y1, y2
pred: 4.0840704496667311000, -0.5877852522924734700, -0.8090169943749467900
real: 4.0840704496667311000, -0.5877852522924732500, -0.8090169943749473400

step 53: t, y1, y2
pred: 4.1626102660064763000, -0.5224985647159489100, -0.8526401643540918500
real: 4.1626102660064763000, -0.5224985647159485800, -0.8526401643540924000

step 54: t, y1, y2
pred: 4.2411500823462216000, -0.4539904997395464700, -0.8910065241883676800
real: 4.2411500823462216000, -0.4539904997395461400, -0.8910065241883682300

step 55: t, y1, y2
pred: 4.3196898986859669000, -0.3826834323650890600, -0.9238795325112867400
real: 4.3196898986859669000, -0.3826834323650886700, -0.9238795325112871800

step 56: t, y1, y2
pred: 4.3982297150257121000, -0.3090169943749462900, -0.9510565162951536400
real: 4.3982297150257121000, -0.3090169943749459000, -0.9510565162951540900

step 57: t, y1, y2
pred: 4.4767695313654574000, -0.2334453638559038600, -0.9723699203976766700
real: 4.4767695313654574000, -0.2334453638559034200, -0.9723699203976771100

step 58: t, y1, y2
pred: 4.5553093477052027000, -0.1564344650402289000, -0.9876883405951377700
real: 4.5553093477052027000, -0.1564344650402284300, -0.9876883405951381000

step 59: t, y1, y2
pred: 4.6338491640449480000, -0.0784590957278425430, -0.9969173337331278500
real: 4.6338491640449480000, -0.0784590957278420290, -0.9969173337331281900

step 60: t, y1, y2
pred: 4.7123889803846932000, 0.0000000000000028311, -0.9999999999999996700
real: 4.7123889803846932000, 0.0000000000000033690, -1.0000000000000000000

step 61: t, y1, y2
pred: 4.7909287967244385000, 0.0784590957278481770, -0.9969173337331274100
real: 4.7909287967244385000, 0.0784590957278487460, -0.9969173337331276300

step 62: t, y1, y2
pred: 4.8694686130641838000, 0.1564344650402344800, -0.9876883405951368800
real: 4.8694686130641838000, 0.1564344650402350600, -0.9876883405951371000

step 63: t, y1, y2
pred: 4.9480084294039290000, 0.2334453638559093600, -0.9723699203976753400
real: 4.9480084294039290000, 0.2334453638559099700, -0.9723699203976755600

step 64: t, y1, y2
pred: 5.0265482457436743000, 0.3090169943749516700, -0.9510565162951518700
real: 5.0265482457436743000, 0.3090169943749522800, -0.9510565162951519800

step 65: t, y1, y2
pred: 5.1050880620834196000, 0.3826834323650942800, -0.9238795325112845200
real: 5.1050880620834196000, 0.3826834323650949400, -0.9238795325112846300

step 66: t, y1, y2
pred: 5.1836278784231649000, 0.4539904997395515200, -0.8910065241883651200
real: 5.1836278784231649000, 0.4539904997395521300, -0.8910065241883651200

step 67: t, y1, y2
pred: 5.2621676947629101000, 0.5224985647159536800, -0.8526401643540888500
real: 5.2621676947629101000, 0.5224985647159543500, -0.8526401643540888500

step 68: t, y1, y2
pred: 5.3407075111026554000, 0.5877852522924780200, -0.8090169943749434500
real: 5.3407075111026554000, 0.5877852522924786900, -0.8090169943749433400

step 69: t, y1, y2
pred: 5.4192473274424007000, 0.6494480483301885400, -0.7604059656000262700
real: 5.4192473274424007000, 0.6494480483301892100, -0.7604059656000261600

step 70: t, y1, y2
pred: 5.4977871437821459000, 0.7071067811865523500, -0.7071067811865421300
real: 5.4977871437821459000, 0.7071067811865530100, -0.7071067811865420200

step 71: t, y1, y2
pred: 5.5763269601218912000, 0.7604059656000355900, -0.6494480483301775500
real: 5.5763269601218912000, 0.7604059656000362600, -0.6494480483301774400

step 72: t, y1, y2
pred: 5.6548667764616365000, 0.8090169943749517800, -0.5877852522924663600
real: 5.6548667764616365000, 0.8090169943749525600, -0.5877852522924661400

step 73: t, y1, y2
pred: 5.7334065928013818000, 0.8526401643540962900, -0.5224985647159413600
real: 5.7334065928013818000, 0.8526401643540969500, -0.5224985647159411400

step 74: t, y1, y2
pred: 5.8119464091411270000, 0.8910065241883715600, -0.4539904997395386400
real: 5.8119464091411270000, 0.8910065241883722300, -0.4539904997395383100

step 75: t, y1, y2
pred: 5.8904862254808723000, 0.9238795325112899600, -0.3826834323650809600
real: 5.8904862254808723000, 0.9238795325112906200, -0.3826834323650805700

step 76: t, y1, y2
pred: 5.9690260418206176000, 0.9510565162951562000, -0.3090169943749380100
real: 5.9690260418206176000, 0.9510565162951567500, -0.3090169943749375100

step 77: t, y1, y2
pred: 6.0475658581603629000, 0.9723699203976785600, -0.2334453638558954000
real: 6.0475658581603629000, 0.9723699203976791100, -0.2334453638558948400

step 78: t, y1, y2
pred: 6.1261056745001081000, 0.9876883405951389900, -0.1564344650402203200
real: 6.1261056745001081000, 0.9876883405951395500, -0.1564344650402197100

step 79: t, y1, y2
pred: 6.2046454908398534000, 0.9969173337331284100, -0.0784590957278338970
real: 6.2046454908398534000, 0.9969173337331288500, -0.0784590957278332310

step 80: t, y1, y2
pred: 6.2831853071795987000, 0.9999999999999995600, 0.0000000000000114908
real: 6.2831853071795987000, 1.0000000000000000000, 0.0000000000000121896

标签: none

评论已关闭