MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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]



  • Prev by Date: Help wanted to find out if bounding function is pierced for n even > 10^7.
  • Next by Date: curved von Koch
  • Previous by thread: Using StoppingTest in NDSolve
  • Next by thread: [Bug] NullSpace[m, Method->OneStepRowReduction] can cause crashes in v5