Re: Using StoppingTest in NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg49620] Re: [mg49533] Using StoppingTest in NDSolve
- From: Ramesh Raju Mudunuri <rmudunuri at uh.edu>
- Date: Sat, 24 Jul 2004 03:48:45 -0400 (EDT)
- References: <200407230959.FAA20289@smc.vnet.net>
- Reply-to: Ramesh Raju Mudunuri <rmudunuri at uh.edu>
- Sender: owner-wri-mathgroup at wolfram.com
Hi,
It looks like that StoppingTest does not take delayed values as
part of the test. It could be due to the way NDSolve stores values
of dependent variables as it integrates.
Since the derivatives (in time) vanish when the dependent variable
or variables attain a constant value, we can use this fact to stop NDSolve
from integrating further-as you wanted. I assume that we are solving an
IVP of order two or more. It appears that StoppingTest only takes current
values of dependent variables as part of the test. We cannot use y'[x]<1
as a test, but can use y[x]>3 as part of the test. If we want to use derivatives
as in y'[x]<1 as a part of the test, then we need to define another function
v(x)=f'(x) and use v[x] in the test. Since an IVP can be written as a set of
first order IVPs in the derivatives of the dependant variables, we can input
the original IVP as a set of first order IVPs to NDSolve. This would allow us
to use derivatives in the StoppingTest.
Coming to your first example f''[x] + f[x] == 0 can be written as set of two
coupled ODEs as {v'[x]+f[x]==0, f'[x]==v[x]}. Instead of solving for f[x], we
now solve for f[x] and v[x] which is f'[x] at each step of integration. StoppingTest
can now be used as
sol = NDSolve[{v'[x] == -y[x], y'[x] == v[x], y[0] == 0, v[0] == 1}, {y,
v}, {x, 0, 10}, StoppingTest -> (Abs[v[x]] < 0.1)][[1]]
{y, v} = {y, v} /. sol
tmax = y[[1]][[1, 2]] (NDSolve integrated from 0 to tmax)
Plot[y[x], {x, 0, tmax}, PlotRange -> All]
Plot[{v[x], y'[x]}, {x, 0, tmax}, PlotRange -> All]
Note that the performance of the StoppingTest depends on the x values
which depend on time steps that NDSolve takes. Look at
http://support.wolfram.com/mathematica/kernel/Symbols/System/StoppingTest.html.
Your problem can be solved as (derx(t)=x'(t) and dery(t)=y'(t))
R = 0.15;
x1 = Sqrt[3] - 1; x2 = -Sqrt[3] - 1; x3 = 0 - 1;
y1 = 1; y2 = 1; y3 = -2;
c = 0.2;
d = 0.25;
magnetx1 = Derivative[1][derx][t] +
R*Derivative[1][x][t] - (x1 - x[t])/
Sqrt[(x1 - x[t])^2 + (y1 - y[t])^2 + d^2]^3 -
(x2 - x[t])/Sqrt[(x2 - x[t])^2 + (y2 - y[t])^2 +
d^2]^3 - (x3 - x[t])/
Sqrt[(x3 - x[t])^2 + (y3 - y[t])^2 + d^2]^3 +
c*x[t] == 0
magnety1 =
Derivative[1][dery][t] + R*Derivative[1][y][t] -
(y1 - y[t])/Sqrt[(x1 - x[t])^2 + (y1 - y[t])^2 +
d^2]^3 - (y2 - y[t])/
Sqrt[(x2 - x[t])^2 + (y2 - y[t])^2 + d^2]^3 -
(y3 - y[t])/Sqrt[(x3 - x[t])^2 + (y3 - y[t])^2 +
d^2]^3 + c*y[t] == 0
solution =
NDSolve[{magnetx1, magnety1, x'[t] == derx[t], y'[t] == dery[t],
x[0] == 2.5,
derx[0] == 0, y[0] == 2.5,
dery[0] == 0}, {x, y, derx, dery}, {t, 0, 100},
MaxSteps -> 10000,
StoppingTest ->
If[t < 11, False, (Abs[derx[t]] < 0.1 && Abs[dery[t]] < 0.1)]][[1]]
OR
solution =
NDSolve[{magnetx1, magnety1, x'[t] == derx[t], y'[t] == dery[t],
x[0] == 2.5,
derx[0] == 0, y[0] == 2.5,
dery[0] == 0}, {x, y, derx, dery}, {t, 0, 100},
MaxSteps -> 10000,
StoppingTest ->
If[t < 11,
False, (Abs[derx[t]] < 0.1 && Abs[dery[t]] < 0.1 &&
Abs[y[t] - y3] < 0.1 || Abs[y[t] - y1] < 0.1 ||
Abs[y[t] - y2] < 0.1)]][[1]]
{x, y, derx, dery} = {x, y, derx, dery} /. solution
tmax = x[[1]][[1, 2]]
Plot[{x[t], y[t]}, {t, 0, tmax}, PlotRange -> All]
Plot[{derx[t], dery[t]}, {t, 0, tmax}, PlotRange -> All]
Hope this helps.
Ramesh
----- Original Message -----
From: Urijah Kaplan
To: mathgroup at smc.vnet.net
To: mathgroup at smc.vnet.net
Sent: Friday, July 23, 2004 4:59 AM
Subject: [mg49620] [mg49533] Using StoppingTest in NDSolve
Hello,
I don't know if what I want to do is possible...briefly, I have a
numerical diff. equation of a chaotic pendulum, that is I trace the path
of a pendulum over three magnets and the pendulum eventually stops over
one of them. In order to save (computational) time, I would like NDsolve
to stop when it's clear that the pendulum is going to stay over one of the
magnets. I also want it to progress for a minimum of time, say t (the
independent variable) going to 100.
Here is a what I tried for a similar example
NDSolve[{f''[x] + f[x] == 0, f[0] == 0, f'[0] == 1}, f, {
x, 0, 10}, StoppingTest -> (Which[x < 5,
False, Abs[f[x] - f[x - 1]] < .1, True])]
here if x is less than 5 I want to continue; and if 5<x<10, and if the
absolute difference between f[x] and f[x-1] is less than a small amount
(say .1) NDSolve will stop.
I'm pretty sure the problem here is that Mtm doesn't know what to do with
f[x-1]. Does anyone have any suggestions? Also, if someone can figure out
how to do this, how would I save the stopping value to a variable? Say
NDSolve stopped at t=20.2, (i.e. the InterpolatingFunction stops at 20.2)
how do I save that number.
Thank you very much for any help.
--Urijah Kaplan
If it helps, here is the actual problem
R = 0.15;
x1 = Sqrt[3] - 1; x2 = -Sqrt[3] - 1; x3 = 0 - 1;
y1 = 1; y2 = 1; y3 = -2;
c = 0.2;
d = 0.25;
magnetx1 = Derivative[2][x][t] +
R*Derivative[1][x][t] - (x1 - x[t])/
Sqrt[(x1 - x[t])^2 + (y1 - y[t])^2 + d^2]^3 -
(x2 - x[t])/Sqrt[(x2 - x[t])^2 + (y2 - y[t])^2 +
d^2]^3 - (x3 - x[t])/
Sqrt[(x3 - x[t])^2 + (y3 - y[t])^2 + d^2]^3 +
c*x[t] == 0; magnety1 =
Derivative[2][y][t] + R*Derivative[1][y][t] -
(y1 - y[t])/Sqrt[(x1 - x[t])^2 + (y1 - y[t])^2 +
d^2]^3 - (y2 - y[t])/
Sqrt[(x2 - x[t])^2 + (y2 - y[t])^2 + d^2]^3 -
(y3 - y[t])/Sqrt[(x3 - x[t])^2 + (y3 - y[t])^2 +
d^2]^3 + c*y[t] == 0;
solution = NDSolve[{magnetx1, magnety1, x[0] == 2.5,
Derivative[1][x][0] == 0, y[0] == 2.5,
Derivative[1][y][0] == 0}, {x, y}, {t, 0, 100},
MaxSteps -> 10000];
xy[t_] = {x[t], y[t]} /. solution;
Plot[Evaluate[{x[t], y[t]} /. solution], {t, 0, 50},
PlotRange -> All]
- References:
- Using StoppingTest in NDSolve
- From: Urijah Kaplan <uak@sas.upenn.edu>
- Using StoppingTest in NDSolve