NDSolve Singularity in solving LLG equation
- To: mathgroup at smc.vnet.net
- Subject: [mg126321] NDSolve Singularity in solving LLG equation
- From: Jiwan Kim <hwoarang.kim at gmail.com>
- Date: Tue, 1 May 2012 05:21:24 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
Hello, Mathgroup.
I am trying to solve the famous Landau-Lifshitz-Gilbert equation, which
explains the dynamic motion of the magnetization in ferromagnetics.
(dM/dt = -(gamma*mu0)MxHeff + alpha/M (MxdM/dt). (x : vector product)
The reducing magnitude of magnetization at time=0 triggers the precession
motion. Using this equation I want to get the theta[t] and phi[t] in
spherical coordinate system.
But, now I have a problem that the error message "NDSolve::ndsz: At t ==
-2.1091991871443846`*^-13, step size is effectively zero; singularity or
stiff system suspected".
I don't think that the equation includes the diverging point. And this
error message always appears when the value of 'SpinTempData' starts to
increase from 300.
I suspect that D[Mr[t], t] is not working well in this code. I don't know
the reason.
Could you help me? Plz.....
In this code, I used (r, theta, phi) parameter instead of (x,y,z), but
Cartesian basis vector was used when the vector product was conducted (in
eq2 expression in the code,
I used 'x' instead of CrossProduct to get a vectorproduct in Cartesian
coordinate)
I attached the following code (SpinTempData is included) :
---------------------------------------------------------------------------------------------------------------------------
Remove["Global`*"];
Needs["VectorAnalysis`"]
SetCoordinates[Spherical];
\[Mu]0 = 4 \[Pi] 10^-7;(* [H/m] *) Ms =
0.66/\[Mu]0;(* 525 emu/cc=0.66 *) Hext =
0.45/\[Mu]0; (* Ms, Hext [A/m] *)
\[Alpha]f = 0.05; \[Gamma] = 1.76 10^11; K0 = -0.5 10^4; \[Psi] =
35 \[Pi]/180.0; ps = 10^-12; \[Beta] = 1.34 10^-5;
\[Nu] = 0.3; B =
1.8 10^11; \[Lambda]s = -3.4 10^-5; Tc = 631; Tp = 98; Rac = -0.39; \
\[Sigma]s1 = 1;
\[Xi]1 = 13.5; \[Xi]2 = 14.5; Nx = 1/3; Ny = 1/3; Nz =
1/3; \[Rho] = 8910; v = 4082;
SpinTempData = {{-5.`, 300}, {-4.95`, 300}, {-4.9`, 300}, {-4.85`,
300}, {-4.8`, 300}, {-4.75`, 300}, {-4.7`, 300}, {-4.65`,
300}, {-4.6`, 300}, {-4.55`, 300}, {-4.5`, 300}, {-4.45`,
300}, {-4.4`, 300}, {-4.35`, 300}, {-4.3`, 300}, {-4.25`,
300}, {-4.2`, 300}, {-4.15`, 300}, {-4.1`, 300}, {-4.05`,
300}, {-4.`, 300}, {-3.95`, 300}, {-3.9`, 300}, {-3.85`,
300}, {-3.8`, 300}, {-3.75`, 300}, {-3.7`, 300}, {-3.65`,
300}, {-3.6`, 300}, {-3.55`, 300}, {-3.5`, 300}, {-3.45`,
300}, {-3.4`, 300}, {-3.35`, 300}, {-3.3`, 300}, {-3.25`,
300}, {-3.2`, 300}, {-3.15`, 300}, {-3.1`, 300}, {-3.05`,
300}, {-3.`, 300}, {-2.95`, 300}, {-2.9`, 300}, {-2.85`,
300}, {-2.8`, 300}, {-2.75`, 300}, {-2.7`, 300}, {-2.65`,
300}, {-2.6`, 300}, {-2.55`, 300}, {-2.5`, 300}, {-2.45`,
300}, {-2.4`, 300}, {-2.35`, 300}, {-2.3`, 300}, {-2.25`,
300}, {-2.2`, 300}, {-2.15`, 300}, {-2.1`, 300}, {-2.05`,
300}, {-2.`, 300}, {-1.95`, 300}, {-1.9`, 300}, {-1.85`,
300}, {-1.8`, 300}, {-1.75`, 300}, {-1.7`, 300}, {-1.65`,
300}, {-1.6`, 300}, {-1.55`, 300}, {-1.5`, 300}, {-1.45`,
300}, {-1.4`, 300}, {-1.35`, 300}, {-1.3`, 300}, {-1.25`,
300}, {-1.2`, 300}, {-1.15`, 300}, {-1.1`, 300}, {-1.05`,
300}, {-1.`, 300}, {-0.95`, 300}, {-0.9`, 300}, {-0.85`,
300}, {-0.8`, 300}, {-0.75`, 300}, {-0.7`, 300}, {-0.65`,
300}, {-0.6`, 300}, {-0.55`, 300}, {-0.5`, 300}, {-0.45`,
300}, {-0.4`, 300}, {-0.35`, 300.00029`}, {-0.3`,
300.00057`}, {-0.25`, 300.00588`}, {-0.2`, 300.05339`}, {-0.15`,
300.48233`}, {-0.1`, 303.32062`}, {-0.05`,
314.09299`}, {-9.39526`*^-15, 337.36538`}, {0.05`,
369.43267`}, {0.1`, 402.14368`}, {0.15`, 430.44968`}, {0.2`,
453.56358`}, {0.25`, 472.37784`}, {0.3`, 487.89433`}, {0.35`,
500.84204`}, {0.4`, 511.7486`}, {0.45`, 521.00431`}, {0.5`,
528.89807`}, {0.55`, 535.64818`}, {0.6`, 541.43276`}, {0.65`,
546.38944`}, {0.7`, 550.63141`}, {0.75`, 554.25236`}, {0.8`,
557.32945`}, {0.85`, 559.92689`}, {0.9`, 562.10125`}, {0.95`,
563.89989`}, {1.`, 565.36016`}, {1.05`, 566.51965`}, {1.1`,
567.4031`}, {1.15`, 568.04153`}, {1.2`, 568.45306`}, {1.25`,
568.66236`}, {1.3`, 568.68462`}, {1.35`, 568.53867`}, {1.4`,
568.23797`}, {1.45`, 567.79709`}, {1.5`, 567.22686`}, {1.55`,
566.53977`}, {1.6`, 565.74438`}, {1.65`, 564.85152`}, {1.7`,
563.86864`}, {1.75`, 562.80459`}, {1.8`, 561.66652`}, {1.85`,
560.46151`}, {1.9`, 559.19632`}, {1.95`, 557.87689`}, {2.`,
556.50931`}, {2.05`, 555.09879`}, {2.1`, 553.6507`}, {2.15`,
552.16966`}, {2.2`, 550.66042`}, {2.25`, 549.1271`}, {2.3`,
547.57386`}, {2.35`, 546.0044`}, {2.4`, 544.42235`}, {2.45`,
542.83101`}, {2.5`, 541.23351`}, {2.55`, 539.6328`}, {2.6`,
538.03158`}, {2.65`, 536.43244`}, {2.7`, 534.83766`}, {2.75`,
533.24951`}, {2.8`, 531.66993`}, {2.85`, 530.10083`}, {2.9`,
528.54384`}, {2.95`, 527.00058`}, {3.`, 525.47245`}, {3.05`,
523.96065`}, {3.1`, 522.4663`}, {3.15`, 520.99046`}, {3.2`,
519.53391`}, {3.25`, 518.09775`}, {3.3`, 516.68226`}, {3.35`,
515.28843`}, {3.4`, 513.91649`}, {3.45`, 512.56702`}, {3.5`,
511.24031`}, {3.55`, 509.93651`}, {3.6`, 508.65601`}, {3.65`,
507.39864`}, {3.7`, 506.16475`}, {3.75`, 504.95404`}, {3.8`,
503.76671`}, {3.85`, 502.60244`}, {3.9`, 501.46125`}, {3.95`,
500.34284`}, {4.`, 499.24705`}, {4.05`, 498.17363`}, {4.1`,
497.12225`}, {4.15`, 496.09268`}, {4.2`, 495.08452`}, {4.25`,
494.09751`}, {4.3`, 493.13119`}, {4.35`, 492.18529`}, {4.4`,
491.25934`}, {4.45`, 490.353`}, {4.5`, 489.46584`}, {4.55`,
488.59747`}, {4.6`, 487.74747`}, {4.65`, 486.91542`}, {4.7`,
486.10092`}, {4.75`, 485.30354`}, {4.8`, 484.52288`}, {4.85`,
483.75851`}, {4.9`, 483.01004`}, {4.95`, 482.27705`}, {5.`,
481.55915`}, {5.5`, 475.12589`}, {6.`, 469.80928`}, {6.5`,
465.32763`}, {7.`, 461.47028`}, {7.5`, 458.08379`}, {8.`,
455.05749`}, {8.5`, 452.31161`}, {9.`, 449.78841`}, {9.5`,
447.44569`}, {10.`, 445.25233`}, {10.5`, 443.18502`}, {11.`,
441.22612`}, {11.5`, 439.36201`}, {12.`, 437.58201`}, {12.5`,
435.87759`}, {13.`, 434.24183`}, {13.5`, 432.669`}, {14.`,
431.15426`}, {14.5`, 429.69348`}, {15.`, 428.28306`}, {15.5`,
426.91983`}, {16.`, 425.60098`}, {16.5`, 424.32397`}, {17.`,
423.08651`}, {17.5`, 421.8865`}, {18.`, 420.72202`}, {18.5`,
419.5913`}, {19.`, 418.49269`}, {19.5`, 417.42466`}, {20.`,
416.38581`}, {20.5`, 415.37479`}, {21.`, 414.39037`}, {21.5`,
413.43139`}, {22.`, 412.49675`}, {22.5`, 411.58543`}, {23.`,
410.69646`}, {23.5`, 409.82893`}, {24.`, 408.98198`}, {24.5`,
408.1548`}, {25.`, 407.34662`}, {25.5`, 406.55672`}, {26.`,
405.78441`}, {26.5`, 405.02904`}, {27.`, 404.28998`}, {27.5`,
403.56666`}, {28.`, 402.85852`}, {28.5`, 402.16502`}, {29.`,
401.48566`}, {29.5`, 400.81996`}, {30.`, 400.16746`}}
datanum = Count[SpinTempData[[All, 1]], _Real];
Curie[temp_] :=
Ms ((1 - Exp[(temp - Tc)/Tc])/(1 - Exp[(300 - Tc)/Tc]))^0.51;
SpinData =
Table[{ps SpinTempData[[i, 1]], Curie[SpinTempData[[i, 2]]]}, {i, 1,
datanum}];
SpinDataitp = Interpolation[SpinData];
Plot[SpinDataitp[t], {t, -5 ps, 30 ps}, PlotRange -> All]
Mr[t_] := If[t > -5 ps, SpinDataitp[t], Ms];
Plot[Mr[t], {t, -5 ps, 30 ps}, PlotRange -> All]
Module[{t, \[Theta]2, \[Phi]2},
Energy[t_, \[Theta]2_, \[Phi]2_] := (1/2 \[Mu]0 Nz Mr[
t]^2) Cos[\[Theta]2]^2 +
1/2 \[Mu]0 Mr[
t]^2 Sin[\[Theta]2]^2 (Ny Sin[\[Phi]2]^2 +
Nx Cos[\[Phi]2]^2) - \[Mu]0 Hext Mr[
t] (Sin[\[Theta]2] Sin[\[Phi]2] Sin[\[Psi]] +
Cos[\[Theta]2] Cos[\[Psi]]);
\[Theta]0 =
FindArgMin[
Energy[-5 ps, \[Theta]2, \[Pi]/2], {\[Theta]2, \[Pi]/4}];];
\[Theta]0 = \[Theta]0[[1]]
M = {Mr[t] Sin[\[Theta][t]] Cos[\[Phi][t]],
Mr[t] Sin[\[Theta][t]] Sin[\[Phi][t]], Mr[t] Cos[\[Theta][t]]};
Heff1 = {(-Mr[t] (Nz Cos[\[Theta][t]]^2 +
Ny Sin[\[Theta][t]]^2 Sin[\[Phi][t]]^2 +
Nx Sin[\[Theta][t]]^2 Cos[\[Phi][t]]^2) +
Hext (Sin[\[Theta][t]] Sin[\[Phi][t]] Sin[\[Psi]] +
Cos[\[Theta][t]] Cos[\[Psi]])) Sin[\[Theta][t]] Cos[\[Phi][
t]] + ((Nz Mr[
t] - (2 K0 + 3 \[Lambda]s \[Sigma]s1)/(\[Mu]0 Mr[t]) -
Mr[t] Ny Sin[\[Phi][t]]^2 -
Mr[t] Nx Cos[\[Phi][t]]^2) Sin[\[Theta][t]] Cos[\[Theta][
t]] + Hext (Cos[\[Theta][t]] Sin[\[Phi][t]] Sin[\[Psi]] -
Sin[\[Theta][t]] Cos[\[Psi]])) Cos[\[Theta][t]] Cos[\[Phi][
t]] - (Mr[
t] (Nx - Ny) Sin[\[Theta][t]] Sin[\[Phi][t]] Cos[\[Phi][t]] +
Hext Cos[\[Phi][t]] Sin[\[Psi]]) Sin[\[Phi][
t]], (-Mr[t] (Nz Cos[\[Theta][t]]^2 +
Ny Sin[\[Theta][t]]^2 Sin[\[Phi][t]]^2 +
Nx Sin[\[Theta][t]]^2 Cos[\[Phi][t]]^2) +
Hext (Sin[\[Theta][t]] Sin[\[Phi][t]] Sin[\[Psi]] +
Cos[\[Theta][t]] Cos[\[Psi]])) Sin[\[Theta][t]] Sin[\[Phi][
t]] + ((Nz Mr[
t] - (2 K0 + 3 \[Lambda]s \[Sigma]s1)/(\[Mu]0 Mr[t]) -
Mr[t] Ny Sin[\[Phi][t]]^2 -
Mr[t] Nx Cos[\[Phi][t]]^2) Sin[\[Theta][t]] Cos[\[Theta][
t]] + Hext (Cos[\[Theta][t]] Sin[\[Phi][t]] Sin[\[Psi]] -
Sin[\[Theta][t]] Cos[\[Psi]])) Cos[\[Theta][t]] Sin[\[Phi][
t]] + (Mr[
t] (Nx - Ny) Sin[\[Theta][t]] Sin[\[Phi][t]] Cos[\[Phi][t]] +
Hext Cos[\[Phi][t]] Sin[\[Psi]]) Cos[\[Phi][t]],
(-Mr[t] (Nz Cos[\[Theta][t]]^2 +
Ny Sin[\[Theta][t]]^2 Sin[\[Phi][t]]^2 +
Nx Sin[\[Theta][t]]^2 Cos[\[Phi][t]]^2) +
Hext (Sin[\[Theta][t]] Sin[\[Phi][t]] Sin[\[Psi]] +
Cos[\[Theta][t]] Cos[\[Psi]])) Cos[\[Theta][
t]] - ((Nz Mr[
t] - (2 K0 + 3 \[Lambda]s \[Sigma]s1)/(\[Mu]0 Mr[t]) -
Mr[t] Ny Sin[\[Phi][t]]^2 -
Mr[t] Nx Cos[\[Phi][t]]^2) Sin[\[Theta][t]] Cos[\[Theta][
t]] + Hext (Cos[\[Theta][t]] Sin[\[Phi][t]] Sin[\[Psi]] -
Sin[\[Theta][t]] Cos[\[Psi]])) Sin[\[Theta][t]]};
eq1 = {D[Mr[t], t] Sin[\[Theta][t]] Cos[\[Phi][t]] +
D[\[Theta][t], t] Mr[t] Cos[\[Theta][t]] Cos[\[Phi][t]] -
D[\[Phi][t], t] Mr[t] Sin[\[Theta][t]] Sin[\[Phi][t]],
D[Mr[t], t] Sin[\[Theta][t]] Sin[\[Phi][t]] +
D[\[Theta][t], t] Mr[t] Cos[\[Theta][t]] Sin[\[Phi][t]] +
D[\[Phi][t], t] Mr[t] Sin[\[Theta][t]] Cos[\[Phi][t]],
D[Mr[t], t] Cos[\[Theta][t]] -
D[\[Theta][t], t] Mr[t] Sin[\[Theta][t]]};
eq2 = -\[Gamma] \[Mu]0 M\[Cross]Heff1 + \[Alpha]f/Mr[t] M\[Cross]eq1;
eq3 = \[Theta][-5 ps] == \[Theta]0;
eq4 = \[Phi][-5 ps] == \[Pi]/2;
eq5 = eq1[[2]] == eq2[[2]];
eq6 = eq1[[3]] == eq2[[3]];
solution1 =
NDSolve[{eq5, eq6, eq3, eq4}, {\[Theta][t], \[Phi][t]}, {t, -5 ps,
30 ps}, MaxSteps -> Infinity, MaxStepSize -> {0.05 ps}]