谁能确认我（欧拉法）对这些常微分方程的实现？

Question

我正在尝试使用 Euler 方法在 C 中实现这些方程（Matsuoka 的振荡器 - 此处有详细信息 http://www.ecila.org/ecila_files/content/papers/ACEICMC05.pdf）。我意识到 Euler 不是最准确的方法，但它在我开发时消除了一些复杂性，一旦我完成了这项工作，我计划切换到 Runge Kutta。这是我的实现，但我无法复制论文中的结果，所以我确定我有问题。但是代码翻了多少遍都找不到。

是否有人能够查看此代码并确认我的工作，或发现任何错误？

#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);

}

Answer 1

宏的定义

#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 可能会被计算两次。在这种情况下应该使用函数。

谁能确认我（欧拉法）对这些常微分方程的实现？

Can anyone confirm my (eulers method) implementation of these ordinary differential equations?

c

math

ode