ODE solution problem
- To: mathgroup at smc.vnet.net
- Subject: [mg93105] ODE solution problem
- From: Virgil Stokes <vs at it.uu.se>
- Date: Sun, 26 Oct 2008 01:29:39 -0500 (EST)
I do not understand why the following code (using Mathematica vers. 5.2)
does not give the correct solution. That is, the variable swt in the
equation module is never set to 1, although the value of y(t) exceeds
the value of yhi. The correct solution should have the following appearance,
*
* * *
* *
*
*
where it should reach a maximum value of 0.1 (yhi) at about 15 seconds,
and then start to decrease, reaching a minimum of 0.025 (ylo) at about
30 seconds, etc. That is, the solution should oscillate between ylo and
yhi. The following code however, when plotted shows a linear increase in
y(t). I would appreciate it greatly if someone can point out my error(s)
in the code, since I have not been able to identify why this code does
not work properly.
ClassicalRungeKutta /:
NDSolve`InitializeMethod[ClassicalRungeKutta, __] :=
ClassicalRungeKutta[];
ClassicalRungeKutta[___]["Step"[f_, t_, h_, y_, yp_]]:=
Block[{deltay, k1, k2, k3, k4},
k1 = yp;
k2 = f[t + 1/2 h, y + 1/2 h k1];
k3 = f[t + 1/2 h, y + 1/2 h k2];
k4 = f[t + h, y + h k3];
deltay = h (1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4);
{h, deltay}
];
c = 0.6; g = 9.81;
R = 0.05; r = 0.007; yhi = 0.1; ylo = 0.025;
Qin = 0.00005;
eq[y_,swtch_] := (Qin - swtch*c*Sqrt[2*g*y]*Pi*r^2)/(Pi*R^2);
Remove[equation];
swt = 0;
equation[y_] := Module[{},
If[y <= ylo,swt=0;ex = eq[y,swt]];
If[y >= yhi,swt=1;ex = eq[y,swt]];
ex];
sol = NDSolve[{y'[t]==
equation[y[t]],y[0]==0},y[t],{t,0,35},Method->ClassicalRungeKutta,StartingStepSize->0.02];
--V. Stokes