Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: ODE solution problem

• To: mathgroup at smc.vnet.net
• Subject: [mg93312] Re: [mg93123] Re: ODE solution problem
• From: Virgil Stokes <vs at it.uu.se>
• Date: Mon, 3 Nov 2008 05:29:00 -0500 (EST)

Here is a follow-up to the "ODE solution problem" that I first posted on 
2008-10-26 [mg93105]. Note, if you have trouble with "seeing" some of 
the characters then try viewing with unicode (UTF-8) character encoding.

First, thanks to Jens-Peer Kuska for some very helpful information (see 
below). I had several errors in the original code that was posted. The 
following Mathematica code works as I intended (at least in vers. 5.2) 
for this problem.

*******************************************

c = 0.6; g = 9.81;
R = 0.05; r = 0.007; yhi = 0.1; ylo = 0.025;
Qin = 0.00005;

eq[y_?NumericQ, swtch_?NumericQ] := (Qin - 
swtch*c*Sqrt[2*g*y]*Pi*r2)/(Pi*R2);
swt = 0;

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

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

tsmax = 100;
tpmax = 100;
h = 0.05;
sol = NDSolve[{y'[t]== 
equation[y[t]],y[0]==0},y[t],{t,0,tsmax},Method->ClassicalRungeKutta,StartingStepSize->h]; 

Plot[y[t] /. sol[[1]], {t, 0, tpmax}];

************************************************
I have a short document that gives some background on this problem and a 
more detailed Mathematica notebook that I would be glad to email to 
anyone interested.

--V. Stokes


Jens-Peer Kuska wrote:
> 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*r2)/(Pi*R2);
>>
>> 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
>>
>

• Follow-Ups:
  • Re: Re: Re: ODE solution problem
    • From: Curtis Osterhoudt <cfo@lanl.gov>