1
$\begingroup$

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!

$\endgroup$
9
  • 1
    $\begingroup$ You might want to look at the implementation of MonitorMethod here. $\endgroup$ – J. M.'s ennui Jun 24 '15 at 3:34
  • 1
    $\begingroup$ 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. $\endgroup$ – mikuszefski Jun 24 '15 at 13:15
  • 3
    $\begingroup$ 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? $\endgroup$ – mikuszefski Jun 24 '15 at 13:18
  • 2
    $\begingroup$ Concerning memorization, this is good to ensure that the value does not vary upon one step $t$. Note that you will also get the same values in consecutive simulations assuming you hit the same $t$'s, which I am not sure about. My question is the following (as I cannot test your code at the moment): Do you see noise on the result? Does it change if you change the temperature (independent of the question if the result is correct with respect to the time step)? $\endgroup$ – mikuszefski Jun 24 '15 at 15:47
  • 1
    $\begingroup$ Well, try to have a look at some old threads, at the documentation center $\endgroup$ – mikuszefski Jun 26 '15 at 11:43
1
$\begingroup$

It seems that the modeling of the step input for the current is rather not physical. The invariant is delayed perturbed in a non-neglectable manner.

First, take into account that NDSolve is proven in practice and not a built-in that is made for step size control. That was a wish by the Mathematica customers and a consequence of the test suits on which it is verified for the highest possible versatility. It saves time and effort for the users.

There are still of course workarounds to control step size for tedious and error-prone partial differential equations. This case is not close to such problems.

I did put the precision options out of the command. Automatic means for the Projection option too the use of Runge-Kutta. A proper choice for it is to match the order of the coefficients favored. So with Automatic these whole lot of optimization are done by Mathematica in NDSolve for the best solution.

With this in mind the results for

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;
sol = 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 -> 1, 
  Method -> {"Projection", Method -> "Automatic", 
    "Invariants" -> (mX[t]^2 + mY[t]^2 + mZ[t]^2 == 1)}(*,
  StepMonitor\[RuleDelayed]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]]]*)]

are

InterpolationFunction

Graphics

The peak at the delayed time t about 5.5 ns for the invariant is a violation of the foundational physics for the problem. Since the problem is reduced from the mathematical standpoint this is the solution.

The best path for a more physical solution will be to model that the peak does not appear. For example model some counter current from the other side of the electric system or some immanent current sources that delay in the calculated or near calculated manner.

Another interpretation might be, that there are times of unphysical solutions unavoidable in this modeling. Another one is a sign is set wrong. Nevertheless, the system recovers.

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;
sol0 = 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]
StepDataPlot[sol0, PlotRange -> {tSt, tFin}]
Plot[{mX[t], mY[t], mZ[t], mX[t]^2 + mY[t]^2 + mZ[t]^2} /. sol0, {t, 
  tSt, tFin}]

Output

Output

Great

$\endgroup$

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.