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>
- Re: Re: ODE solution problem