       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},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},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