我试图以 Octave 为单位求解一个微分方程,但我选择的微分单位需要很长时间,所以我决定用 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
我选择了 I get a segmentation fault;是什么原因造成的?
备用 Octave 解决方案
在我的一个 friend 告诉我他能够使用 MATLAB 使用 1e-3
的微分步求解方程后,我发现 MATLAB 有一个“硬”版本的 ode23
模块——“stiff”的意思是专门用于求解那些需要极小步长的微分方程。后来我在 Octave 中搜索了“刚性”ODE 求解器,发现 lsode
属于该类别。第一次尝试时,lsode
以微秒的速度求解了方程(均比 MATLAB 和我的 C 实现更快),并且得到了完美的解。 FOSS 万岁!
最佳答案
你的递归终止速度不够快,所以你在浪费你的堆栈。
要解决这个问题,只需将其设为循环即可,看起来您实际上并没有在做任何需要递归的事情。
我认为这样做:
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)的情况下编译它,任何体面的编译器实际上都足够聪明,可以进行尾递归调用,并且您的程序不会出现段错误。
关于c - 数值微分方程求解器算法意外出现段错误,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16612757/