
Problems with NDSolve
Posted:
Jul 25, 1996 10:39 AM


Mathgroupers,
I have a system of ODE, dx/dt = f(t,x,y), dy/dt = g(t,x,y), x(0) = x0, y(0) = y0, for which NDSolve is giving unexpected results. Moreover, I have run it with no problems with a fixed stepsize RungeKutta package so I think the math is correct.
The functions f and g are piecewise defined (they are both essentially continuous), and it is around where the pieces fit together that NDSolve has problems. If you have some expertise with NDSolve and/or variable stepsize algorithms could you take a look at the system and see if you can spot what is happening?
The system:
f(t,x,y) = If[ Sqrt[9x^4+12x^3+4x^2+4y^2] >= 0.1, Sign[x]*2y/Sqrt[9x^4+12x^3+4x^2+4y^2], 1/Sqrt[3x+2] ]
g(t,x,y) = If[ Sqrt[9x^4+12x^3+4x^2+4y^2] >= 0.1, Sign[x]*(3x^2+2x)/Sqrt[9x^4+12x^3+4x^2+4y^2], Sqrt[3x+1]/Sqrt[3x+2] ]
x0 = 1, y0 = 0
The trajectory of this lies on the algebraic curve y^2 = x^2 + x^3. This curve has a double point at (0,0) which is where the second definitions in f and g kick in. To keep this simple I have given you the system so that smoothness of solution is guaranteed only for t >= 0.
The problem with NDSolve is that it "bounces back" from the origin at (0,0), instead of passing through it as the fixed step RungeKutta method does. What follows is the NDSolve version with all the variables printed out so that you can see it "bouncing back" at the origin. (Sorry about the length of this, but I have found that NDSolve is picky about how you input a piecewise defined equation. I need the flexibility of using a Module in order to obtain the desired solution for all t.)
f[t_,x_,y_,defn_]:= Module[{a}, Print["Using piecewise definition ",defn,". Variable t = ",t]; Print["(x,y) = ",{x,y}]; If[ defn==1, a=Sign[x]*2y/Sqrt[9x^4+12x^3+4x^2+4y^2], a=1/Sqrt[3x+2] ]; Print["dx/dt = ",a]; Return[a]; ];
g[t_,x_,y_,defn_]:= Module[{a}, If[ defn==1, a=Sign[x]*(3x^2+2x)/Sqrt[9x^4+12x^3+4x^2+4y^2], a=Sqrt[3x+1]/Sqrt[3x+2] ]; Print["dy/dt = ",a];Print[""]; Return[a]; ];
dx[t_,x_,y_]:= If[ Sqrt[9x^4+12x^3+4x^2+4y^2 ] < 0.1, f[t,x,y,2], f[t,x,y,1] ];
dy[t_,x_,y_]:= If[ Sqrt[9x^4+12x^3+4x^2+4y^2 ] < 0.1, g[t,x,y,2], g[t,x,y,1] ];
NDSolve[ {x'[t] == dx[t,x[t],y[t]], y'[t] == dy[t,x[t],y[t]], x[0] == 1, y[0]== 0}, {x,y},{t,0,2}]
Thanks for any help you can give me.
Steve Wilkinson wilkinson@nku.edu

