Re: Modification of Variable in NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg107923] Re: Modification of Variable in NDSolve
- From: dh <dh at metrohm.com>
- Date: Wed, 3 Mar 2010 05:49:28 -0500 (EST)
- References: <hmiikh$3r6$1@smc.vnet.net>
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 at metrohm.com>
Internet:<http://www.metrohm.com>