1
$\begingroup$

I want to solve a second-order ODE in the form of y^{''} = \frac{a (y^{'})^2}{b y^{'}+cy+d} by numerical method (eg, solver ODE45), given initial condition of $y(0)$ and y'(0). The results are wield and numbers go out of machinery bound. I guess the catch is that what is in the denominator becomes highly unstable when it converges to zero. I tried bound it away from zero with no avail.

Could anyone provide insights on how to proceed with the numerical procedure? Thanks in advance...

  • 0
    Fine, `R`. B$u$t no code?2011-11-29

1 Answers 1

1

Discretizing the ODE by finite differences gives

$\frac{y_2-2y_1 + y_0}{h^2} = \frac{a\left(\frac{y_1-y_0}{h}\right)^2}{b\left(\frac{y_1-y_0}{h}\right) + cy_1 + d},$ or $y_2 = 2y_1 - y_0 + \frac{ah(y_1-y_0)^2}{b(y_1-y_0)+chy_1+dh}.$ Here's C++ code I wrote which has no trouble integrating this ODE for $a=-10,b=c=d=1$, initial conditions y(0)=0,y'(0)=10 and time step $h=0.01$. I'm sure you can adapt it to whatever language you prefer:

#include   using namespace std;  double y0, y1;  void step(double dt) {     double y2 = 2*y1-y0 - 10*dt*(y1-y0)*(y1-y0)/(y1-y0+dt*y1+dt);     y0 = y1;     y1 = y2; }  int main() {     y0 = 0;     y1 = .1;     for(int i=0; i<1000; i++)     {         step(0.01);         cout << y1 << endl;     } }