Prev: Efficient representation of many different primitives
Next: Harmonic Numbers
From: Shawn Garbett on 2 Mar 2010 03:35
I need to modify a variable while a system of ODE's are being solved.
I've simplified the problem to the following:

dA = B[t];
dB = (1 - A[t]^2) B[t] - A[t];
dF = B[t] + 0.2;

Module[{count = 0, last = 0.},
sol = NDSolve[{A'[t] == dA, B'[t] == dB, F'[t] == dF, A[0] == 1,
B[0] == 1, F[0] == 1}, {A, B, F}, {t, 0, 100},
StepMonitor :> If[last <= 5 <= F[t] && last < F[t],
count++; last = F[t],
last = F[t]]
]; count]

This does a threshold crossing of 5 for the variable F from below. The
count that results is correct, but one more thing is needed. This
shows the graphs

disp = Plot[#[t] /. sol, {t, 0, 100}, PlotRange -> All] &;
GraphicsGrid[Map[disp, {{A, B}, {F}}, {2}]]

What I need is for F to be set to half it's value at each threshold
crossing from below. I've been trying introducing another variable,
but so far nothing I've tried seems to work.

From: Patrick Scheibe on 2 Mar 2010 08:00
Hi,

I didn't read you post carefully but it seems you want to catch an event during
the run of NDSolve and maybe this helps you out.

http://reference.wolfram.com/mathematica/tutorial/NDSolveEventLocator.html

Cheers
Patrick


Am Mar 2, 2010 um 9:35 AM schrieb Shawn Garbett:

> I need to modify a variable while a system of ODE's are being solved.
> I've simplified the problem to the following:
>
> dA = B[t];
> dB = (1 - A[t]^2) B[t] - A[t];
> dF = B[t] + 0.2;
>
> Module[{count = 0, last = 0.},
> sol = NDSolve[{A'[t] == dA, B'[t] == dB, F'[t] == dF, A[0] == 1,
> B[0] == 1, F[0] == 1}, {A, B, F}, {t, 0, 100},
> StepMonitor :> If[last <= 5 <= F[t] && last < F[t],
> count++; last = F[t],
> last = F[t]]
> ]; count]
>
> This does a threshold crossing of 5 for the variable F from below. The
> count that results is correct, but one more thing is needed. This
> shows the graphs
>
> disp = Plot[#[t] /. sol, {t, 0, 100}, PlotRange -> All] &;
> GraphicsGrid[Map[disp, {{A, B}, {F}}, {2}]]
>
> What I need is for F to be set to half it's value at each threshold
> crossing from below. I've been trying introducing another variable,
> but so far nothing I've tried seems to work.
>


From: dh on 3 Mar 2010 05:49
Hi,
if you Change the function midway it will no more be continuous, not to
speak of differentiable. Therefore, you need to restart the Solver at
each discontinuity.
Here is a somewhat laborious method:
A and B are independent of F. Therefore, first solve for A and B.
F only depends on B:
dA = B[t];
dB = (1 - A[t]^2) B[t] - A[t];
dF = B[t] + 0.2;

funB = B /.
NDSolve[{A'[t] == dA, B'[t] == dB, A[0] == 1, B[0] == 1}, {A,
B}, {t, 0, 100}][[1]]
For F we then have a linear first order DE:
F'[t]==0.2+B[t]
Because of linearity you can solve:
F'[t]==0.2
F'[t]=B[t]
separately:
f1 = F /. NDSolve[{F'[t] == 0.2, F[0] == 0}, F, {t, 0, 100}][[1]]
f2 = F /. NDSolve[{F'[t] == funB[t], F[0] == 1}, F, {t, 0, 100}][[1]]
The first ODE results in a linear growth, the second in an oscillation.
Our function is therefore:
fun1[t_]:=f1[t]+f2[t]
Now where do we hit 5 from below for the first time? Near 20
Plot[{fun1[t]-5},{t,0,100}]
or numerically:
r1=t/. FindRoot[fun1[t] - 5, {t, 20}]
the searched function becomes now:
fun2[t_] := fun1[t] - If[t > r1, 0.5, 0] fun1[r1]
Now where does this hit 5 from below, near 33:
Plot[fun2[t] - 5, {t, 0, 100}]
r2=t/. FindRoot[fun2[t] - 5, {t, 33}]
fun3[t_] := fun2[t] - If[t > r2, 0.5, 0] fun2[r1]
We continue:
Plot[fun3[t] - 5, {t, 0, 100}]
r3=t/. FindRoot[fun3[t] - 5, {t, 40}]
fun4[t_] := fun3[t] - If[t > r2, 0.5, 0] fun3[r1]

I hope you got the idea and you can continue.
Daniel

On 02.03.2010 09:35, Shawn Garbett wrote:
> I need to modify a variable while a system of ODE's are being solved.
> I've simplified the problem to the following:
>
> dA = B[t];
> dB = (1 - A[t]^2) B[t] - A[t];
> dF = B[t] + 0.2;
>
> Module[{count = 0, last = 0.},
> sol = NDSolve[{A'[t] == dA, B'[t] == dB, F'[t] == dF, A[0] == 1,
> B[0] == 1, F[0] == 1}, {A, B, F}, {t, 0, 100},
> StepMonitor :> If[last<= 5<= F[t]&& last< F[t],
> count++; last = F[t],
> last = F[t]]
> ]; count]
>
> This does a threshold crossing of 5 for the variable F from below. The
> count that results is correct, but one more thing is needed. This
> shows the graphs
>
> disp = Plot[#[t] /. sol, {t, 0, 100}, PlotRange -> All]&;
> GraphicsGrid[Map[disp, {{A, B}, {F}}, {2}]]
>
> What I need is for F to be set to half it's value at each threshold
> crossing from below. I've been trying introducing another variable,
> but so far nothing I've tried seems to work.
>


--

Daniel Huber
Metrohm Ltd.
Oberdorfstr. 68
CH-9100 Herisau
Tel. +41 71 353 8585, Fax +41 71 353 8907
E-Mail:<mailto:dh(a)metrohm.com>
Internet:<http://www.metrohm.com>


From: Shawn Garbett on 5 Mar 2010 04:32
I took these suggestion into account and now have the following, it's
crude but works. If I can refine it into a loop, then I'll be done.

First Interval

x1 = Reap[NDSolve[
{A'[t] == dA, B'[t] == dB, F'[t] == dF,
A[0] == 1, B[0] == 1, F[0] == 1},
{A, B, F}, {t, 0, 100},
Method -> {
"EventLocator",
"Event" -> F[t] - 5,
"Direction" -> 1,
"EventAction" :> Sow[t]
}
]]

sol1 = x1[[1]][[1]]

div1 = x1[[2]][[1]][[1]]

Second Interval

x2 = Reap[NDSolve[
{A'[t] == dA, B'[t] == dB, F'[t] == dF,
A[div1] == (A[div1] /. sol1),
B[div1] == (B[div1] /. sol1),
F[div1] == (F[div1] /. sol1)/2},
{A, B, F}, {t, div1, 100},
Method -> {
"EventLocator",
"Event" -> F[t] - 5,
"Direction" -> 1,
"EventAction" :> Sow[t]
}
]]

sol2 = x2[[1]][[1]]

div2 = x2[[2]][[1]][[1]]

Third Interval

x3 = Reap[NDSolve[
{A'[t] == dA, B'[t] == dB, F'[t] == dF,
A[div2] == (A[div2] /. sol2),
B[div2] == (B[div2] /. sol2),
F[div2] == (F[div2] /. sol2)/2},
{A, B, F}, {t, div2, 100},
Method -> {
"EventLocator",
"Event" -> F[t] - 5,
"Direction" -> 1,
"EventAction" :> Sow[t]
}
]]

sol3 = x3[[1]][[1]]

div3 = x3[[2]][[1]][[1]]

AA[t_ /; t <= div1] := (A[t] /. sol1);
AA[t_ /; div1 < t <= div2] := (A[t] /. sol2) ;
AA[t_ /; div2 < t <= div3] := (A[t] /. sol3);

BB[t_ /; t <= div1] := (B[t] /. sol1);
BB[t_ /; div1 < t <= div2] := (B[t] /. sol2) ;
BB[t_ /; div2 < t <= div3] := (B[t] /. sol3);

FF[t_ /; t <= div1] := (F[t] /. sol1);
FF[t_ /; div1 < t <= div2] := (F[t] /. sol2) ;
FF[t_ /; div2 < t <= div3] := (F[t] /. sol3);

disp = Plot[#[t], {t, 0, div3}, PlotRange -> All] &;
GraphicsGrid[Map[disp, {{AA, BB}, {FF}}, {2}]]


From: DrMajorBob on 7 Mar 2010 05:10
Looks like a lot of great code, there, but it doen't work, not a bit...
since dA, dB, and dF are undefined.

Bobby

On Fri, 05 Mar 2010 03:31:55 -0600, Shawn Garbett
<shawn.garbett(a)gmail.com> wrote:

> I took these suggestion into account and now have the following, it's
> crude but works. If I can refine it into a loop, then I'll be done.
>
> First Interval
>
> x1 = Reap[NDSolve[
> {A'[t] == dA, B'[t] == dB, F'[t] == dF,
> A[0] == 1, B[0] == 1, F[0] == 1},
> {A, B, F}, {t, 0, 100},
> Method -> {
> "EventLocator",
> "Event" -> F[t] - 5,
> "Direction" -> 1,
> "EventAction" :> Sow[t]
> }
> ]]
>
> sol1 = x1[[1]][[1]]
>
> div1 = x1[[2]][[1]][[1]]
>
> Second Interval
>
> x2 = Reap[NDSolve[
> {A'[t] == dA, B'[t] == dB, F'[t] == dF,
> A[div1] == (A[div1] /. sol1),
> B[div1] == (B[div1] /. sol1),
> F[div1] == (F[div1] /. sol1)/2},
> {A, B, F}, {t, div1, 100},
> Method -> {
> "EventLocator",
> "Event" -> F[t] - 5,
> "Direction" -> 1,
> "EventAction" :> Sow[t]
> }
> ]]
>
> sol2 = x2[[1]][[1]]
>
> div2 = x2[[2]][[1]][[1]]
>
> Third Interval
>
> x3 = Reap[NDSolve[
> {A'[t] == dA, B'[t] == dB, F'[t] == dF,
> A[div2] == (A[div2] /. sol2),
> B[div2] == (B[div2] /. sol2),
> F[div2] == (F[div2] /. sol2)/2},
> {A, B, F}, {t, div2, 100},
> Method -> {
> "EventLocator",
> "Event" -> F[t] - 5,
> "Direction" -> 1,
> "EventAction" :> Sow[t]
> }
> ]]
>
> sol3 = x3[[1]][[1]]
>
> div3 = x3[[2]][[1]][[1]]
>
> AA[t_ /; t <= div1] := (A[t] /. sol1);
> AA[t_ /; div1 < t <= div2] := (A[t] /. sol2) ;
> AA[t_ /; div2 < t <= div3] := (A[t] /. sol3);
>
> BB[t_ /; t <= div1] := (B[t] /. sol1);
> BB[t_ /; div1 < t <= div2] := (B[t] /. sol2) ;
> BB[t_ /; div2 < t <= div3] := (B[t] /. sol3);
>
> FF[t_ /; t <= div1] := (F[t] /. sol1);
> FF[t_ /; div1 < t <= div2] := (F[t] /. sol2) ;
> FF[t_ /; div2 < t <= div3] := (F[t] /. sol3);
>
> disp = Plot[#[t], {t, 0, div3}, PlotRange -> All] &;
> GraphicsGrid[Map[disp, {{AA, BB}, {FF}}, {2}]]
>
>


--
DrMajorBob(a)yahoo.com

 |  Next  |  Last
Pages: 1 2
Prev: Efficient representation of many different primitives
Next: Harmonic Numbers