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}}]
Improve this question
$\endgroup$
4
  • $\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, 2015 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, 2015 at 16:21
  • $\begingroup$ You know Gamma is a Function, see Gamma, and your [ ] are possibly ill placed. $\endgroup$
    – user9660
    Nov 6, 2015 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, 2015 at 18:31

2 Answers 2

Reset to default
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 *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ FWIW NDSolve does not diverge like that for large t. (perhaps a bifurcation?) $\endgroup$
    – george2079
    Nov 6, 2015 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, 2015 at 19:06
  • $\begingroup$ @george2079 You are correct. Answer updated. $\endgroup$
    – Jack LaVigne
    Nov 6, 2015 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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.