首页 文章

数值微分方程求解器算法意外地断块

提问于
浏览 次
1

我试图用八度音程求解微分方程，但是我选择的微分单元需要永远，所以我决定用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 我选择了我得到了一个分段错误;什么可能导致这个？

备用Octave解决方案

在我的一个朋友告诉我他能够使用MATLAB的差分步骤解决方程后，我发现MATLAB有一个"stiff"版本的 ode23 模块--"stiff"意味着特殊的解决那些需要非常极端的微分方程小步长 . 后来我在Octave中搜索了"stiff" ODE解算器，发现 lsode 属于该类别 . 在第一次尝试时， lsode 以微秒（比MATLAB和我的C实现更快）解决方程式，并提供了完美的解决方案 . 万岁FOSS！

c math segmentation-fault numerical-methods differential-equations

5 回答

  • 0

    你的递归并没有足够快地终止，所以你正在吹嘘你的筹码 .

    为了解决这个问题，只需使它成为一个循环，看起来你实际上并没有做任何需要递归的事情 .

    我认为这样做：

    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），任何体面的编译器实际上都足够聪明，可以进行尾递归调用，以及你的程序不会是段错误 .

    回复于
  • 1

    那么，你已经对你的实际问题有了很多答案 . 我只想提请你注意另一个lib . 使用boost.odeint，如果你需要一个快速且易于使用的颂歌解算器，它基本上就是最先进的库 . 忘记GSL，Matlab等，Odeint优于所有这些 .

    您的程序将如下所示：

    #include <boost/numeric/odeint.hpp>
using namespace boost::numeric::odeint;

typedef boost::array<double,3> State;

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

void my_ode( State const &s , State &dsdt , double t ) {
   double const &z = s[0],    // this is just a name
                &w = s[1],    // forwarding for better
                &i = s[2];    // readability of the ode
   dsdt[0] = w;
   dsdt[1] = 1./J * ( Ka*i - ko*z - bo*w );
   dsdt[2] = 1./L * ( E - R*i - Kb*w );
}

void printer( State const &s , double t ) {
   std::cout << s[0] << " " << s[1] << " " << s[2] << std::endl;
}

int main() {
   State s = {{ 0, 0, 0 }};
   integrate_const(
      euler<State>() , my_ode , s , 0. , 2. , 1e-6 , printer
   );
}
    回复于
  • 3

    求解器函数递归调用自身 . 有多深？可能是几百万次 . 你的堆栈空间不足吗？每个递归都需要四个双倍空间和一个堆栈帧，这很快就会增加 .

    我建议你将解算器重写为迭代而不是递归函数 .

    回复于
  • 3

    正如@hexist所说，这里不需要递归 .

    堆栈溢出：

    ==4734== Memcheck, a memory error detector
==4734== Copyright (C) 2002-2011, and GNU GPL'd, by Julian Seward et al.
==4734== Using Valgrind-3.7.0 and LibVEX; rerun with -h for copyright info
==4734== Command: ./demo
==4734== 
==4734== Stack overflow in thread 1: can't grow stack to 0x7fe801ff8
==4734== 
==4734== Process terminating with default action of signal 11 (SIGSEGV)
==4734==  Access not within mapped region at address 0x7FE801FF8
==4734==    at 0x40054E: solver (demo.c:18)
==4734==  If you believe this happened as a result of a stack
==4734==  overflow in your program's main thread (unlikely but
==4734==  possible), you can try to increase the size of the
==4734==  main thread stack using the --main-stacksize= flag.
==4734==  The main thread stack size used in this run was 8388608.
==4734== Stack overflow in thread 1: can't grow stack to 0x7fe801fe8
==4734== 
==4734== Process terminating with default action of signal 11 (SIGSEGV)
==4734==  Access not within mapped region at address 0x7FE801FE8
==4734==    at 0x4A226E0: _vgnU_freeres (vg_preloaded.c:58)
==4734==  If you believe this happened as a result of a stack
==4734==  overflow in your program's main thread (unlikely but
==4734==  possible), you can try to increase the size of the
==4734==  main thread stack using the --main-stacksize= flag.
==4734==  The main thread stack size used in this run was 8388608.
==4734== 
==4734== HEAP SUMMARY:
==4734==     in use at exit: 0 bytes in 0 blocks
==4734==   total heap usage: 0 allocs, 0 frees, 0 bytes allocated
==4734== 
==4734== All heap blocks were freed -- no leaks are possible
==4734== 
==4734== For counts of detected and suppressed errors, rerun with: -v
==4734== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 4 from 4)
    回复于
  • 2

    你正在解算器上进行1e6递归调用 . 我猜你用完了 . 尝试在求解器内部循环，更新变量而不是调用函数 .

    伪代码：

    while t < tf
      do dt step
      t = t + dt

    等等

    回复于

相关问题