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>

**Help wanted to find out if bounding function is pierced for n even > 10^7.**

**curved von Koch**

**Using StoppingTest in NDSolve**

**[Bug] NullSpace[m, Method->OneStepRowReduction] can cause crashes in v5**