1
$\begingroup$

I'm trying to find My[t], Mx[t] and Mz[t] values when evaluated at t = Infinity.

Code to date:

γ = 2.675*10^8;
T1 = 100*10^-3;
T2 = 10*10^-3;
B0 = 2;
M0 = 1*10^-1;
B1 = 1*10^-7;

(*Above are constants*)

NDSolve[{Mx'[t] == -Mx[t]/T2, Mx[0] == 0}, Mx, {t, 0, 0.5}]

NDSolve[{My'[t] == γ*Mz[t]*B1 - My[t]/T2,
  Mz'[t] == -γ*My[t]*B1 + (M0 - Mz[t])/T1,
  My[0] == 0, Mz[0] == M0}, {My, Mz}, {t, 0, 0.5}, 
 MaxSteps -> Infinity]

Plot[Evaluate[Mx[t] /. %], {t, 0, 0.5}, 
 FrameLabel -> {Style[t, 14], Style[Subscript[M, x], 12]}, 
 PlotTheme -> "Scientific"]

Plot[Evaluate[{Mz[t], My[t]} /. %], {t, 0, 0.5}, 
 FrameLabel -> {Style[t, 14], Style[M, 14]}, 
 PlotTheme -> "Scientific", 
 PlotLegends -> 
  Placed[{"\!\(\*SubscriptBox[\(M\), \(y\)]\)", 
    "\!\(\*SubscriptBox[\(M\), \(z\)]\)"}, {0.9, 0.85}], 
 ClippingStyle -> None, PlotRange -> All, 
 FrameTicks -> {{Automatic, All}, {Automatic, None}}]
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – user9660 Nov 6 '15 at 16:20
  • $\begingroup$ The styling of your question can be improved. Please take the time to read the markdown help page. You'll better communicate your problem when you use the right formatting. $\endgroup$ – user9660 Nov 6 '15 at 16:21
  • $\begingroup$ You know Gamma is a Function, see Gamma, and your [ ] are possibly ill placed. $\endgroup$ – user9660 Nov 6 '15 at 16:26
  • $\begingroup$ If the code is evaluated in order you are using % incorrectly. (first usage tries to get MX from the My,Mz solution for example) I'd suggest actually assiging the NDSolve results to symbols and using those instead of relying on output history (%) $\endgroup$ – george2079 Nov 6 '15 at 18:31
2
$\begingroup$

Update

george2079 pointed out that using approximate numbers created a problem with the solution for high values of t. Replacing the approximate value for gamma with a rational number, the problem is stable at high values, including infinity.

γ = 267500000;
t1 = 1/10;
t2 = 1/100;
b0 = 2;
m0 = 1/10;
b1 = 1/10000000;

DSolve works with your equations.

The first differential equation yields a constant value for mx[t] of zero so I will skip it.

The second differential equation gives:

solYZ = 
 DSolve[{my'[t] == γ*mz[t]*b1 - my[t]/t2, 
   mz'[t] == -γ*my[t]*b1 + (m0 - mz[t])/t1, my[0] == 0, 
   mz[0] == m0}, {my[t], mz[t]}, t]

(* {{my[t] -> (1/11501679980)
   107 E^(-55 t - (Sqrt[20951] t)/4 - 
     1/4 (-220 + Sqrt[20951]) t) (1676080 E^((Sqrt[20951] t)/2) - 
      838040 E^(1/4 (-220 + Sqrt[20951]) t) - 
      18649 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) - 
      838040 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) + 
      18649 Sqrt[20951]
        E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t)), 
  mz[t] -> (1/11501679980)
   E^(-55 t - (Sqrt[20951] t)/4 - 
     1/4 (-220 + Sqrt[20951]) t) (670432000 E^((Sqrt[20951] t)/2) + 
      239867999 E^(1/4 (-220 + Sqrt[20951]) t) - 
      2518780 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) + 
      239867999 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) + 
      2518780 Sqrt[20951]
        E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t))}} *)

We can create functions out of the solutions in a number of ways. Below is one method:

my[t_] = solYZ[[1, 1, 2]]

(* (1/11501679980)107 E^(-55 t - (Sqrt[20951] t)/4 - 
  1/4 (-220 + Sqrt[20951]) t) (1676080 E^((Sqrt[20951] t)/2) - 
   838040 E^(1/4 (-220 + Sqrt[20951]) t) - 
   18649 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) - 
   838040 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) + 
   18649 Sqrt[20951]
     E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t)) *)

and for mz[t]

mz[t_] = solYZ[[1, 2, 2]]

(* (1/11501679980)E^(-55 t - (Sqrt[20951] t)/4 - 
  1/4 (-220 + Sqrt[20951]) t) (670432000 E^((Sqrt[20951] t)/2) + 
   239867999 E^(1/4 (-220 + Sqrt[20951]) t) - 
   2518780 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) + 
   239867999 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) + 
   2518780 Sqrt[20951]
     E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t)) *)

Now we can plot it:

Plot[{my[t], mz[t]}, {t, 0, 1}, 
 PlotStyle -> {{Thick, Blue}, {Thick, Red}}]

Mathematica graphics

It appears to quickly approach a constant. If you take the limits you will find:

Limit[my[t], t -> ∞]
(* 428/27449 *)

Limit[mz[t], t -> ∞]
(* 1600/27449 *)
$\endgroup$
  • $\begingroup$ FWIW NDSolve does not diverge like that for large t. (perhaps a bifurcation?) $\endgroup$ – george2079 Nov 6 '15 at 18:45
  • $\begingroup$ I'm pretty sure there is a numerical precision issue with evaluation the analytic solution. If you make all the parameters exact ( gamma=267500000) this goes away. $\endgroup$ – george2079 Nov 6 '15 at 19:06
  • $\begingroup$ @george2079 You are correct. Answer updated. $\endgroup$ – Jack LaVigne Nov 6 '15 at 20:10
1
$\begingroup$

Not sure if this is what you are looking for but it's a start:

γ = 2.675*10^8; 
T1 = 100*10^-3; 
T2 = 10*10^-3; 
B0 = 2; 
M0 = 1*10^-1; 
B1 = 1*10^-7;

nds1 = NDSolve[{Mx'[t] == -Mx[t]/T2, Mx[0] == 0}, Mx, {t, 0, 0.5}]

nds2 = NDSolve[{My'[t] == γ*Mz[t]*B1 - My[t]/T2,
Mz'[t] == -γ*My[t]*B1 + (M0 - Mz[t])/T1, My[0] == 0, 
Mz[0] == M0}, {My, Mz}, {t, 0, 0.5}, MaxSteps -> Infinity]

Plot[Evaluate[Mx[t] /. nds1], {t, 0, 0.5}, 
FrameLabel -> {Style[t, 14], Style[Subscript[M, x], 12]}, 
PlotTheme -> "Scientific"]

Plot[Evaluate[{Mz[t], My[t]} /. nds2], {t, 0, 0.5}, 
PlotTheme -> "Scientific", PlotRange -> All, 
FrameTicks -> {{Automatic, All}, {Automatic, None}}]

enter image description here

enter image description here

enter image description here

$\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.