NDSolve breaking down

I'm trying to model a situation involving charged sphere in a dynamic electric potential, and find out how the rotational motion of the sphere affects the translational dynamics in two dimensions. Here q0 refers to charge at the centre of the sphere, whereas q1 is charge deposited at its surface, which will lead to a torque. ωd is the driving frequency of the electric field. My model seems to work well in certain limits, but breaks down in the region I'm actually interested in. My equations are

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
eqns1 = {
x''[t] == - (1/m) E0 Cos[ωd t] (q0 x[t] + q1 Abs[Sin[ϕ[t]]] (x[t] + R Sin[ϕ[t]]))
- γt x'[t],
y''[t] == -(1/m) E0 Cos[ωd t] (q0 y[t] + q1 Abs[Cos[ϕ[t]]] (y[t] + R Cos[ϕ[t]]))
- γt y'[t] - 9.8 m,
ϕ''[t] == 1/MI E0 Cos[ωd t] q1 R (-Cos[ϕ[t]] (x[t] + R Sin[ϕ[t]]) +
Sin[ϕ[t]] (y[t] + R Cos[ϕ[t]])) - γr ϕ'[t],
x'[0] == xv0,
ϕ'[0] == ϕv0,
x[0] == x0,
ϕ[0] == ϕ0,
y'[0] == yv0,
y[0] == x0};


I then have a list of parameters (which I'll at the end of the post), including the initial values. I solve the equations using

sol[p_, time_] :=
NDSolve[eqns1 //. p, {x[t], ϕ[t], ϕ'[t], x'[t], y[t], y'[t]}, {t, 0, time},
MaxSteps -> ∞, "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}];


Where I have experimented with a number of Methods, including "StiffnessSwitching", "BDH" , "ExplicitRungeKutta" and "ImplicitRungeKutta". This does work with certain parameters, basically when q0 is large and q1 is low, but breaks down if I try to switch those values, making q1>q0. I run this with time=1. Here are the parameters which do work:

p = {
m -> 2000 4/3 π R^2,
R -> 10^-7,
E0 -> 10^14,
MI -> 5/2 m R^2,
q0 -> q0b eV ,
q0b -> 100,
q1 -> q1b eV,
q1b -> 5,
eV -> 1.602 10^-19,
γt -> 1,
γr -> 1,
ωd -> 15000,
xv0 -> 10^-2,
ϕv0 -> 1,
x0 -> 10^-9,
y0 -> 10^-9,
yv0 -> 10^-2,
ϕ0 -> π/10};


However it breaks down when I send q0b-> 5, q1b-> 100. The error code I get is:

At t == 0.099279970214417, step size is effectively zero;
singularity or stiff system suspected


Any help would be very much appreciated!

• What is the value of time? Also, Cos[ω d t] should be Cos[ωd t] Nov 16, 2015 at 22:47
• Ok, updated to address these points. Typically running it with time=1. Nov 17, 2015 at 14:30
• Look at your equations numerically (eqns1 //. p // N) and work on normalizing things to get rid of the numerically large/small values. ( introduce tprime=t wd , xprime=X/R etc. ). It may not help, the system may simply be too unstable to integrate. Nov 17, 2015 at 19:57
• @DanGoldwater It turns out that your system of equations is chaotic, and this probably explains the problems we both had solving them. I am sending you this comment, because I prepared an answer, deleted it, edited it, and finally undeleted it. I have learned the hard way that you will not receive any notice of the new version of the answer. Best wishes. Nov 18, 2015 at 2:39
• I think you'd better explain the model a little more, for example, the meaning of all the parameters and the original equation etc.… Well, actually what I want to say is, are you sure the equation set is correct? I'm not familiar with electromagnetism, but I guess the -9.8 is gravity acceleration? If so, then why it's multiplied by m? Nov 18, 2015 at 4:23

In the course of trying to solve this problem, I have run a number of physical and numerical parameters. The most striking behavior that I observed is extreme sensitivity of the results to those parameters. This suggests that the system is chaotic, especially for large q1b. To illustrate this extreme sensitivity, compute ϕ for fixed values of x and y:

sc = NDSolveValue[{ϕ''[t] == 1/MI E0 Cos[ωd t] q1 R (-Cos[ϕ[t]] (x + R Sin[ϕ[t]]) +
Sin[ϕ[t]] (y + R Cos[ϕ[t]])) - γr ϕ'[t], ϕ'[0] == 1, ϕ[0] == π/10} //.
p /. {x -> 4*10^-5, y -> 2*10^-5}, {ϕ, ϕ'}, {t, 0, 1/100}];
ParametricPlot[{sc[[1]][t], sc[[2]][t]/ωd /. p}, {t, 0, 1/100},
AxesLabel -> {ϕ, ϕ'}, AspectRatio -> 1]


with q0b -> 5, q1b -> 100.

In contrast, increasing y to 2.01*10^-5 yields

Thus, although

sol[p_, time_] := NDSolve[eqns1 //. p, {x[t], y[t], φ[t],}, {t, 0, time}, MaxSteps -> ∞,
MaxStepFraction -> 10^-3/(ωd time) /. p, Method -> "StiffnessSwitching"];


gives a stable answer, at least out to time = 1/10`,

changing anything produces a very different answer. Sometimes that different answer is unstable, sometimes not.