# Decrease calculation time of NDSolve for coupled differential equations

Hi Stackexchange community,

I am rather new to Mathematica, especially when it comes to dealing with numerical stuff. Right now I am trying to numerically solve a system of two coupled differential equations for certain sets of values. The system is the following:

    modp =
{D[r[t]*Cos[ϕ[t]], {t, 2}] == w1*D[r[t]*Sin[ϕ[t]], t] - w2*Cos[w1*t]*Sin[ϕ[t]],
D[r[t]*Sin[ϕ[t]], {t, 2}] == -w1*D[r[t]*Cos[ϕ[t]], t] - w2*Cos[w1*t]*Cos[ϕ[t]]};


For the values

v0 = 7.5*10^7;
w1 = 1.52*10^10;
w2 = 1.59*10^(-11);


Mathematica is able to solve these differential equations numerically. But for the values

v0 = 7.5*10^7;
w1 = 1.52*10^10;
w2 = 0.437;


I get the message

NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 2.879453677024192*^-7.

Now, I know that the maximum number of steps can be increased by the MaxSteps-command but if I increase the steps by a factor of 10, the same warning occurs for a time that is a factor of 10 higher than the previous one:

dglmodphotp = NDSolve[
{modp, r == v0/w1, ϕ == 0, r' == 0, ϕ' == -v0/r},
{r[t], ϕ[t]}, {t, 0, Pi*10^7}, MaxSteps -> 100000]


Maximum number of 100000 steps reached at the point t == 2.843461540218415*^-6.

So in order to solve the System of ODE for this time intervall I have to increase the maximum number of steps by a factor of 10^14 which (how I see it) roughly also increases the calculation time by this factor. Is there another possibility to solve this system which is not that time-comsuming ?

EDIT: The value of v0 is 7.5*10^7.

• Hi, welcome to the site. Can you give the value of v0 you are using? Your very large w1 means that the $sin$ and $cos$ terms are oscillating wildly with a period of $10^{-10}$, and then you are trying to integrate up to $10^7$. On my system it claims to integrate, but the results don't actually satisfy the ODEs. – KraZug Jul 9 '18 at 7:33
• Same observations as @KraZug. Maybe the w1 should be w1 = 1.52*10^(-10)  ? – Akku14 Jul 9 '18 at 7:47
• Hi, I edited the value of v0 to be 7.5*10^7. The value of w1 is correct. The ODE system describes the Motion of an electron due to the Lorentz force law plus a small modification. The frequency w1 is e*B/m with the elemental charge e, the magnetic field B and the electron mass m. For a magnetic field of B=0.08 Tesla (it has to have this value) i get this large value of w1. – Levigeddon Jul 9 '18 at 7:54
• Then it should be no surprise that you need a mammoth number of steps to integrate to such a large value of $t$ - the underlying functions are oscillating with such a short timestep and then you are trying to go to huge time. More worrying to me is that the solutions I'm getting out are not consistent with the input equations. – KraZug Jul 9 '18 at 8:13
• To follow up on what @KraZug has mentioned, could you rescale things so that your numbers are more reasonable and less wildly oscillatory? As long as you did it in some consistent, intelligent way I expect you could back the target result out. – b3m2a1 Jul 9 '18 at 8:26

What is expected to get as a result of this task? The function $\phi (t)$ is practically linear, the function $r(t)$ is practically constant.

modp = {D[r[t]*Cos[\[Phi][t]], {t, 2}] ==
w1*D[r[t]*Sin[\[Phi][t]], t] - w2*Cos[w1*t]*Sin[\[Phi][t]],
D[r[t]*Sin[\[Phi][t]], {t, 2}] == -w1*D[r[t]*Cos[\[Phi][t]], t] -
w2*Cos[w1*t]*Cos[\[Phi][t]]};
w1 = 1.52*10^10;
w2 = 0.437; v0 = 7.5*10^7;
sol = NDSolveValue[{modp, r == v0/w1, \[Phi] == 0,
r' == 0, \[Phi]' == -v0/r}, {r[t], \[Phi][t],
r'[t], \[Phi]'[t]}, {t, 0, Pi*10^4}];
{Plot[sol[], {t, 0, Pi*10^4}],
Plot[sol[], {t, 0, Pi*10^4}]} • Yes, that is basically what I expect. r(t) describes the radius of the motion of an electron induced by a magnetic field and a small modification so I don't expect it to Change by a large ammount. Can you tell me the main difference between NDSolve and NDSolveValue ? – Levigeddon Jul 9 '18 at 9:57
• Using NDSolveValue, we get a solution in the form of an interpolation function, so we can use it directly, like any function. – Alex Trounev Jul 9 '18 at 10:02
• Nice solution but you used only 1/1000 of the timerange (OP) – Ulrich Neumann Jul 9 '18 at 10:06
• For the original range I get warnings again but by increasing the MaxSteps und decreasing the MaxStepSize they vanish. Seems like the Change of r(t) is really negligible even for such huge time Intervalls. – Levigeddon Jul 9 '18 at 11:02
• Yes, I got a solution up to $t=\pi *10^4$, but I did not see anything significant - almost a constant and a linear function. What can you expect more? – Alex Trounev Jul 9 '18 at 11:04