# Speed up NDSolve for a changing equation?

For physics research, there is a differential equation that I simulate again and again. It would be wonderful to speed it up. Each time I run it, part of the input function $$I(t)$$ changes. It takes a few minutes each time, and after running a few hundred iterations per day, it is eating up a bunch of time. Each time I run it, the input function $$I(t)$$ only changes a little bit, and the output $$\textbf{m(t)}$$ also only changes a little bit.

The equation is the Landau -Lifshitz-Gilbert-Slonczewski Equation (LLGS) :

$$\frac{d\textbf{m}}{dt}=-|\gamma|\textbf{m} \times \textbf{B}_\textrm{eff}-|\gamma|\alpha_\textrm{g}\textbf{m}\times[\textbf{m}\times\textbf{B}_\textrm{eff}]+|\gamma|\alpha_\textrm{j}~I(t)~\textbf{m}\times[\textbf{m}\times\textbf{p}]$$

In this equation, I am solving for $$\bf{m}=\bf{m}(t)$$, which is magnetization (a 3-D vector). $$I(t)$$ is current, an input parameter that is a function of time, which changes each time the LLGS is solved. $$t$$ for time, $$\gamma$$, $$\alpha_\textrm{g}$$, $$\alpha_\textrm{j}$$ are scalar constants, and $$\textbf{B}_\textrm{eff}$$ and $$\textbf{p}$$ are constant vectors.

Here is my code for solving this equation:

gamma = 176;
alphag = 0.01;
alphajConstant = 0.00603;
p = {Cos[Pi/6], 0, Sin[Pi/6]};
current[t_] := Piecewise[{{3, t <= 50}, {(5-3)/(150-50) (t - 50) + 3, 50 < t < 150}, {5, t >= 150}}] + .01*Sin[2*Pi*30*t];
Beff[t_] := {0, 0, 1.5 - 0.8*(m[t].{0, 0, 1})};
cons[t_] := -gamma*Cross[m[t], Beff[t]];
tGilbdamp[t_] := alphag*Cross[m[t], cons[t]];
tSlondamp[t_] := current[t]*alphajConstant*gamma*Cross[m[t], Cross[m[t], p]];
LLGS = {m'[t] == cons[t] + tGilbdamp[t] + tSlondamp[t], m[0] == {0, 0, 1}};
sol1 = NDSolve[LLGS, {m}, {t, 0, 200}, MaxSteps -> \[Infinity]];
mm[t_] = m[t] /. sol1[[1]] ;


To reiterate, this equation is solving for $$\textbf{m(t)}$$=m[t] with a chaging input parameter $$I(t)$$=current[t] . What might change? Perhaps a second current & simulation would look like this:

current[t_] := Piecewise[{{3.1, t <= 60}, {(5.1-3.1)/(200-60) (t - 60) + 3.1, 60< t < 200}, {5.1, t >= 200}}] + .05*Sin[2*Pi*31*t];
LLGS = {m'[t] == cons[t] + tGilbdamp[t] + tSlondamp[t], m[0] == {0, 0, 1}};
sol1 = NDSolve[LLGS, {m}, {t, 0, 200}, MaxSteps -> \[Infinity]];


What I am considering:

One thought is to use the functionality for NDSolve Reinitialize . Will this require me to change how my LLGS equation is formatted? I assume that to use this I will need to change how current[t] appears in my equation, and have it instead be a set of initial conditions. Is this the best option to speed things up?

Perhaps I should define a method? Perhaps there a way to use previous solutions as a starting point?

Is there a better way to make this go faster? Both individually, and in repetition?

• Please show us a minimal working example, currently there're several parameters missing in your code. Jan 20, 2017 at 11:10
• @xzczd I have added working code. Jan 20, 2017 at 11:42
• I do not think equation processing is a significant issue here. Your system is highly oscillatory and simply requires a huge number of increments. I don't see anything you can do but to brute force crunch out each solution. Jan 20, 2017 at 20:06
• NDSolveReinitialize[] will save significant time only if the time for NDSolveProcessEquations[] is significantly long. I suspect that the time taken for NDSolveIterate[] will be at least a few orders of magnitude greater. In other words, I don't expect reinitialization to save much time. Jan 23, 2017 at 23:56
• Just noticed you've removed the introduction for LLGS from your question, may I ask why? Jan 25, 2017 at 3:41

You can speed up the evaluation a good bit by breaking the solution domain so you get rid of the Piecewise:

gamma = 176;
alphag = 0.01;
alphajConstant = 0.00603;
p = {Cos[Pi/6], 0, Sin[Pi/6]};
current[t_] := 3 + .01*Sin[2*Pi*30*t];
Beff[t_] := {0, 0, 1.5 - 0.8*(m[t].{0, 0, 1})};
cons[t_] := -gamma*Cross[m[t], Beff[t]];
tGilbdamp[t_] := alphag*Cross[m[t], cons[t]];
tSlondamp[t_] :=
current[t]*alphajConstant*gamma*Cross[m[t], Cross[m[t], p]];
LLGS = {m'[t] == cons[t] + tGilbdamp[t] + tSlondamp[t],
m[0] == {0, 0, 1}};
sol1 = NDSolve[LLGS, {m}, {t, 0, 50}, MaxSteps -> \[Infinity]];
current[t_] := (5 - 3)/(150 - 50) (t - 50) + 3 + .01*Sin[2*Pi*30*t];
m50 = (m /. First@sol1)[50];
LLGS = {m'[t] == cons[t] + tGilbdamp[t] + tSlondamp[t], m[50] == m50};
sol2 = NDSolve[LLGS, {m}, {t, 50, 150}, MaxSteps -> \[Infinity]];
current[t_] := 5 + .01*Sin[2*Pi*30*t];
m150 = (m /. First@sol2)[150];
LLGS = {m'[t] == cons[t] + tGilbdamp[t] + tSlondamp[t],
m[150] == m150};
sol3 = NDSolve[LLGS, {m}, {t, 150, 200}, MaxSteps -> \[Infinity]];


This runs in about 3-4 seconds..

plot the third component using ListPlot to plot the actual solution points.

Show[{
ListPlot[
Transpose[{Flatten@(m /. sol1[[1]])["Grid"], (m /. sol1[[1]])[
"ValuesOnGrid"][[All, 3]]}]],
ListPlot[
Transpose[{Flatten@(m /. sol2[[1]])["Grid"], (m /. sol2[[1]])[
"ValuesOnGrid"][[All, 3]]}]],
ListPlot[
Transpose[{Flatten@(m /. sol3[[1]])["Grid"], (m /. sol3[[1]])[
"ValuesOnGrid"][[All, 3]]}]]}, PlotRange -> All]


components 1 and 2..

• This solution works well, and is how I implemented it. Jan 24, 2017 at 21:04

You can also convert the Piecewise[] into terms of UnitStep[], using SimplifyPWToUnitStep, and get a significant speed-up:

current[t_] :=
SimplifyPWToUnitStep@
Piecewise[{
{3, t <= 50},
{.02 (t - 50) + 3, 50 < t < 150},
{5, t >= 150}}] + .01*Sin[2*Pi*30*t];
LLGS = {m'[t] == cons[t] + tGilbdamp[t] + tSlondamp[t],  m[0] == {0, 0, 1}};

sol1 = NDSolve[LLGS, {m}, {t, 0, 200}, MaxSteps -> ∞]; // AbsoluteTiming

(*  {2.24653, Null}  *)
`
• This solution works well. Jan 24, 2017 at 21:04