Numerical differential equation solver algorithm segfaults unexpectedly

Question

I was trying to solve a differential equation in octave but it takes forever with the differential unit I have chosen, so I decided to code it in C. Here is the algorithm:

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

The differential unit being defined as h (as you may see), has to be that small, otherwise the values will get out of hand and the solution will not be correct. Now, with numerically greater values of h, the program goes from start to end with no errors (except for nan values), but with the h I have chosen I get a segmentation fault; what could be causing this?

Alternate Octave solution

After a friend of mine told me he was able to solve the equation using a differential step of 1e-3 using MATLAB, I found out that MATLAB has a "stiff" version of its ode23 module--"stiff" meaning special for solving those differential equations that require an extremely small step size. I later searched for "stiff" ODE solvers in Octave and found that lsode pertains in that category. On the first try, lsode solved the equation in microseconds (both faster than MATLAB and my C implementation), and with a perfect solution. Long live FOSS!

Answer 1

You're recursion isn't terminating fast enough, so you're blowing your stack.

To get around this, just make it a loop, it doesn't look like you're actually doing anything that needs recursion.

I think this does it:

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");
}

---

As a side note, your function is eligible for tail recursion optimizations, so if you compile it with some optimizations turned on (-O2 for example), any decent compiler will actually be smart enough to make that a tail recursive call, and your program will not segfault.

Answer 2

Well, you got already plenty of answers on your actual problem. I just want to draw your attention to another lib. Use boost.odeint, that's basically the state-of-the-art library if you need a fast and easy to use ode solver. Forget about GSL, Matlab, etc., Odeint outperforms all of them.

Your program would then look like this:

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

Answer 3

The solver function calls itself recursively. How deep? Probably several million times. Are you running out of stack space? Each recursion needs space for four doubles and a stack frame, that quickly adds up.

I suggest you rewrite the solver as an iterative instead of recursive function.

Answer 4

As @hexist says, you don't need recursion here.

Stack overflow:

==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)

Answer 5

You are doing 1e6 recursive calls on solver. I guess you run out of stack. Try a loop inside solver with updating your variables instead of recalling the function.

pseudo code:

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

and so on