c - 数值微分方程求解器算法段异常

我试着用倍频程来解一个微分方程，但我所选择的微分单位要花很长时间，所以我决定用C来编码。下面是算法：

#include <stdio.h>

double J = 5.78e-5; // (N.m)/(rad/s^2)
double bo = 6.75e-4; // (N.m)/(rad/s)
double ko = 5.95e-4; // (N.m)/rad
double Ka = 1.45e-3; // (N.m)/A
double Kb = 1.69e-3; // V/(rad/s)
double L = 0.311e-3; // mH
double R = 150; // ohms
double E = 5; // V

// Simulacion
int tf = 2;
double h = 1e-6;

double dzdt, dwdt, didt;

void solver(double t, double z, double w, double i) {
    printf("%f %f %f\n", z, w, i);
    if (t >= tf) {
        printf("Finished!\n");
        return; // End simulation
    }
    else {
        dzdt = w;
        dwdt = 1/J*( Ka*i - ko*z - bo*w );
        didt = 1/L*( E - R*i - Kb*w );
        // Solve next step with newly calculated "initial conditions"
        solver(t+h, z+h*dzdt, w+h*dwdt, i+h*didt);
    }
}

int main() {
    solver(0, 0, 0, 0);
    // Solve data
    // Write data to file
    return 0;
}

被定义为h的微分单元（如您所见）必须很小，否则数值将失控，解将不正确现在，当数值大于h时，程序从头到尾都没有错误（除了nan值），但是当我选择h时，我得到了一个分段错误；这是什么原因造成的？
交替倍频程解
在我的一个朋友告诉我他可以使用MATLAB用一个1e-3的微分步来解方程之后，我发现MATLAB有一个ode23模块的“stiff”版本--“stiff”意思是专门用来解那些需要非常小步长的微分方程。我后来在八度音阶中搜索“stiff”的ODE解算器，发现lsode属于这一类在第一次尝试中，lsode以微秒的速度（比MATLAB和我的C实现都快）解出了方程，并且得到了一个完美的解福斯万岁！

最佳答案

你的递归还没有足够快的结束，所以你是吹你的堆栈。
为了解决这个问题，只需将它设为一个循环，它看起来不像是在做任何需要递归的事情。
我想是这样的：

void solver(double t, double z, double w, double i) {
    while (!(t >= tf)) {
        printf("%f %f %f\n", z, w, i);
        dzdt = w;
        dwdt = 1/J*( Ka*i - ko*z - bo*w );
        didt = 1/L*( E - R*i - Kb*w );
        // Solve next step with newly calculated "initial conditions"
        t = t+h;
        z = z+h*dzdt;
        w = w+h*dwdt;
        i = i+h*didt;
    }
    printf("Finished!\n");
}

顺便说一下，您的函数有资格进行尾部递归优化，因此如果您在启用某些优化的情况下编译它（例如（-O2），任何一个好的编译器实际上都足够聪明，能够进行尾部递归调用，并且您的程序不会出现segfult。

关于c - 数值微分方程求解器算法段异常，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/16612757/