MathGroup Archive 2008

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

Search the Archive

Re: ODE solution problem

  • To: mathgroup at smc.vnet.net
  • Subject: [mg93123] Re: ODE solution problem
  • From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
  • Date: Mon, 27 Oct 2008 03:13:37 -0500 (EST)
  • References: <ge0t7t$785$1@smc.vnet.net>

Hi,

how do you think should Mathematica find out:

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];

that y[t] is smaller or larger than ylo/yhi
y[t] is a symbolic expression as to see from here


sol = NDSolve[{y'[t]==
equation[y[t]],y[0]==0},y[t],{t,0,35},Method->ClassicalRungeKutta,StartingStepSize->0.02];


Do you mean

equation[y_?NumericQ] := Module[{},
    If[y <= ylo,swt=0;ex = eq[y,swt]];
    If[y >= yhi,swt=1;ex = eq[y,swt]];
ex];

?? and for what is Module[{},..] good for ?
and is ex not always eq[y,swt] and
is

eq[y,Which[y<ylow,swt=0,y>yhi,swt=1,True,swt]]

not shorter?

Regards
   Jens


Virgil Stokes wrote:
> 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: nested derivatives
  • Next by Date: Re: Re: Function pure for Select
  • Previous by thread: ODE solution problem
  • Next by thread: My talk "The Joy of Tagging" at this year's International