NDSolve breaking down

Question

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["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
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!

Answer 1

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.