Problems with NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg4504] Problems with NDSolve
- From: wilkinson at nku.edu (Steven Wilkinson)
- Date: Fri, 2 Aug 1996 02:22:33 -0400
- Sender: owner-wri-mathgroup at wolfram.com
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 Runge-Kutta 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 Runge-Kutta 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 at nku.edu
==== [MESSAGE SEPARATOR] ====