# Retrieving forthcoming step size from NDSolve

I'm solving a system of ODEs with a random component, which changes at each calculation step and depends on the step size. So I'm trying to realize something like this: (1) size of the forthcoming step is retrieved from NDSolve (2) random component is calculated (3) only then the step is evaluated with a corresponding random component. Unfortunately, I don't know how to squeeze between the moment, when NDSolve selects step size and when it evaluates the step. StepMonitor (commented in the code) doesn't allow that, as I understand, since it returns information after the step is evaluated.

ps1 Fixed step method is 5x longer, which means hours of difference, since I need 10^5-10^6 evaluations for a single final plot, so this would be a really undesirable option.

ps2 Using previous step size also doesn't seem good due to abrupt changes at some points

Here is the code:

Remove["Global*"]

(* PARAMETERS *)
γ =
1.76*10.^7*10.^4 (*T^-1s^-1*);    α = 0.06;   θsh = \
0.3;   Ms = 7.*10.^5; (*A/m*)   tFM = 2.*10.^-9; (*m*)
currAMP = 30.*10.^-3 (*A*);   sHM =
10.*10.^-6*3.*10.^-9(* m^2, tHM = 3 nm *);
Hamp = 0.03; Hk = 0.1;
a = 0. °; b = 0. °; c = 0. °; Ang =
0. °;
temp = 3; rad = 20.*10.^-9; (*Thermal effect*)
kTL = 2.;
tSt = 0.;  tFin = 1.5*10.0^-8;  tStep = 1.*10.0^-12;

(* CURRENT PULSE *)
tPulse = 5.*10^-9; tRise = 0.5*10^-9;
jcurr[t_] := Piecewise[{{1., -1 <= t <= tPulse}}, 0.];
Cjconst = (3.29375*^-16) (θsh/Ms/tFM)*currAMP/sHM;
Cj[t_] := Cj[t] = jcurr[t]*Cjconst;

(* FIELD COMPONENTS *)
σf[temp_, rad_] := Sqrt[(
tStep)];  σ = σf[temp, rad];  (*stand.dev*)
RandComp[] :=
RandomVariate[
NormalDistribution[0., σ]]/Sqrt[3.]; (*random component*)
Happ = ({
{Cos[c], -Sin[c], 0.},
{Sin[c], Cos[c], 0.},
{0., 0., 1.}
}).({
{Cos[b], 0., Sin[b]},
{0., 1., 0.},
{-Sin[b], 0., Cos[b]}
}).({
{1., 0., 0.},
{0., Cos[a], -Sin[a]},
{0., Sin[a], Cos[a]}
}).( {
{Hamp*Cos[Ang]},
{Hamp*(-Sin[Ang])},
{0.}
}); (*permanent applied component*)
Hrand[t_] := Hrand[t] = {RandComp[], RandComp[], RandComp[]};
H[t_, mz_] := Flatten[Happ + Hrand[t] + {0., -kTL*Cj[t], Hk*mz}]

(* ODE SYSTEM *)
t0 = 0;
NDSolve[( {
{mX'[t] == γ (-mY[t]*H[t, mZ[t]][[3]] +
mZ[t]*H[t, mZ[t]][[
2]] - α (mY[
t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]]) -
mZ[t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]])) +
Cj[t]*mY[t]*mX[t] + α*Cj[t]*mZ[t]), mX[0] == 0.},
{mY'[t] == γ (-mZ[t]*H[t, mZ[t]][[1]] +
mX[t]*H[t, mZ[t]][[
3]] - α (mZ[
t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]]) -
mX[t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]])) -
Cj[t] (mX[t]^2 + mZ[t]^2)), mY[0] == 0.},
{mZ'[t] == γ (-mX[t]*H[t, mZ[t]][[2]] +
mY[t]*H[t, mZ[t]][[
1]] - α (mX[
t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]]) -
mY[t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]])) +
Cj[t]*mY[t]*mZ[t] - α*Cj[t]*mX[t]), mZ[0] == 1.}
} ), {mX, mY, mZ}, {t, tSt, tFin}, StartingStepSize -> 10.^-16,
Method -> {"Projection", Method -> "Automatic",
"Invariants" -> (mX[t]^2 + mY[t]^2 + mZ[t]^2 == 1)},
AccuracyGoal -> 15, PrecisionGoal -> 10000 (*,StepMonitor:>Block[{},
δt=t-t0;t0=t;
Print["Current time: ",t];
Print["Current step: ",δt];
Print["mX: ",mX[t]];
Print["mY: ",mY[t]];
Print["mZ: ",mZ[t],"
"];
AbsoluteTiming[Pause[1]]
]*)
]


Thanks!

• You might want to look at the implementation of MonitorMethod here. – J. M. will be back soon Jun 24 '15 at 3:34
• Despite the fact that this is an interesting question, there is a reason why most people use a fixed time step Heun scheme in LLG. So before wasting your time, have a look at this and the references therein. – mikuszefski Jun 24 '15 at 13:15
• Moreover, in a variable time step method you need to decide which step to take. This is usually done by comparing one step with $t$ to two steps with $t/2$. As you have a stochastic component you cannot compare the two endpoints. From this point of view I am surprised that this works in Mathematica. Are you sure that the stochastic field is really generated for every step? – mikuszefski Jun 24 '15 at 13:18
• @mikuszefski, so, to determine standard deviation, size step is required and to select the size step, stochastic component (i.e. standard deviation) should be fixed. So it's a vicious circle. Thank you for that insight and for the reference. Seems like I have to stick to fixed step method. – Alexander Jun 24 '15 at 15:21
• @mikuszefski As for the stochastic field generation, I use memoization for that, so it should be calculated once and then kept for all calculations on that step, until the time t is changed. At least, multiple trials' results are different, so I believe it works. – Alexander Jun 24 '15 at 15:27