10
$\begingroup$

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?

$\endgroup$
7
  • 3
    $\begingroup$ Please show us a minimal working example, currently there're several parameters missing in your code. $\endgroup$
    – xzczd
    Jan 20, 2017 at 11:10
  • $\begingroup$ @xzczd I have added working code. $\endgroup$
    – axsvl77
    Jan 20, 2017 at 11:42
  • $\begingroup$ 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. $\endgroup$
    – george2079
    Jan 20, 2017 at 20:06
  • $\begingroup$ NDSolve`Reinitialize[] will save significant time only if the time for NDSolve`ProcessEquations[] is significantly long. I suspect that the time taken for NDSolve`Iterate[] will be at least a few orders of magnitude greater. In other words, I don't expect reinitialization to save much time. $\endgroup$
    – Michael E2
    Jan 23, 2017 at 23:56
  • $\begingroup$ Just noticed you've removed the introduction for LLGS from your question, may I ask why? $\endgroup$
    – xzczd
    Jan 25, 2017 at 3:41

2 Answers 2

11
$\begingroup$

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]

enter image description here

components 1 and 2..

enter image description here

$\endgroup$
1
  • $\begingroup$ This solution works well, and is how I implemented it. $\endgroup$
    – axsvl77
    Jan 24, 2017 at 21:04
13
$\begingroup$

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

current[t_] := 
  Simplify`PWToUnitStep@
    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}  *)
$\endgroup$
1
  • $\begingroup$ This solution works well. $\endgroup$
    – axsvl77
    Jan 24, 2017 at 21:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.