Re: Re: Re: ODE solution problem

*To*: mathgroup at smc.vnet.net*Subject*: [mg93330] Re: [mg93312] Re: [mg93123] Re: ODE solution problem*From*: Curtis Osterhoudt <cfo at lanl.gov>*Date*: Tue, 4 Nov 2008 06:16:50 -0500 (EST)*Organization*: LANL*References*: <200811031029.FAA05248@smc.vnet.net>*Reply-to*: cfo at lanl.gov

Hi, Virgil, I'd be interested---if you're amenable---in getting a copy of the = notebook(s) you mention in your message. Best wishes, Curtis O. On Monday 03 November 2008 03:29:00 Virgil Stokes wrote: > 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 =E2=89=A4 ylo, swt = 0]; > If[y =E2=89=A5 yhi, swt = 1]; > ex = eq[y, swt]]; > > ClassicalRungeKutta /: > NDSolve`InitializeMethod[ClassicalRungeKutta, __] := ClassicalRungeKutt= a[]; > > 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,S= tartingStepSize->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 > >> > > > > > =2D- ========================== ========================== ======== Curtis Osterhoudt cfo at remove_this.lanl.and_this.gov PGP Key ID: 0x4DCA2A10 Please avoid sending me Word or PowerPoint attachments See http://www.gnu.org/philosophy/no-word-attachments.html ========================== ========================== ========

**References**:**Re: Re: ODE solution problem***From:*Virgil Stokes <vs@it.uu.se>