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[(
2.*α*1.38*10.^-23*temp)/(γ*Ms*tFM*3.1415*rad^2*
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!
MonitorMethodhere. $\endgroup$ – J. M.'s ennui♦ Jun 24 '15 at 3:34