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
    So, MATLAB. It'd be nice if you'd mentioned the MATLAB code you already have.2011-11-29
  • 0
    Actually no. I used R's solver. Basically the function is defined as $$ dg = y[2]; dg^{2} = \frac{a(y[2])^2}{b*y[2]+c*y[1]+d}. $$2011-11-29
  • 0
    Fine, `R`. But 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;
    }
}