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
>