c - 谁能证实我对这些常微分方程的(欧拉方法)实现？

我试图用欧拉方法在C语言中实现这些方程（松冈振子-这里有详细信息）。我意识到欧拉不是最精确的方法，但它在我开发的时候消除了一些复杂性，一旦我有了这个工作，我就打算切换到龙格库塔。
这是我的实现，但我不能复制论文中的结果，所以我肯定我有问题。但不管我翻了多少遍代码，还是找不到。
是否有人可以查看此代码并确认我的工作，或发现任何错误？

#define POSPART(X)  X > 0.0 ? X : 0.0

double matsuoka_calc_nextVal(double in, double t1, double t2,
                          double c, double b, double g,
                          double *x1, double *x2,
                          double *v1, double *v2, double step)
{
    double posX1 = POSPART(*x1);
    double posX2 = POSPART(*x2);
    double posIn = POSPART(in);


    // calculate derivatives
    double dx1 = (c - *x1 - (b*(*v1)) - (g*posX2) - posIn) / t1;
    double dx2 = (c - *x2 - (b*(*v2)) - (g*posX1) - posIn) / t1;
    double dv1 = (posX1 - *v1) / t2;
    double dv2 = (posX2 - *v2) / t2;

    // increment value by 1 euler step (using eulerStep =0.2 in testing)
    *x1 += dx1 * eulerStep;
    *x2 += dx2 * eulerStep;
    *v1 += dv1 * eulerStep;
    *v2 += dv2 * eulerStep;
    return POSPART(*x1) - POSPART(*x2);

}

最佳答案

宏的定义

#define POSPART(X)  X > 0.0 ? X : 0.0

是错误的。
例如，如果你写

POSPART(1.0) - POSPART(1.0)

扩展为

1.0 > 0.0 ? 1.0 : 0.0 - 1.0 > 0.0 ? 1.0 : 0.0

这相当于

(1.0 > 0.0) ? 1.0 : (((0.0 - 1.0) > 0.0) ? 1.0 : 0.0)

因此，POSPART(1.0) - POSPART(1.0)被评估为1.0，因为1.0 > 0.0是真的。
若要避免此类问题，应将参数和函数样式宏的整个表达式括起来，如下所示：

#define POSPART(X)  ((X) > 0.0 ? (X) : 0.0)

在这种情况下，这应该有效，但是像POSPART(x += 1.0)这样的东西不会很好地工作，因为x += 1.0可能会被评估两次。
在这种情况下应该使用函数。