MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Two questions about DSolve
  • Next by Date: Re: Export ListPlots using PlotMarkers to .pdf
  • Previous by thread: Re: Two questions about DSolve
  • Next by thread: Re: ODE solution problem