Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2012

[Date Index] [Thread Index] [Author Index]

Search the Archive

Tracking Differential Equations

  • To: mathgroup at smc.vnet.net
  • Subject: [mg128232] Tracking Differential Equations
  • From: William Duhe <williamduhe at hotmail.com>
  • Date: Thu, 27 Sep 2012 22:47:47 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • Delivered-to: l-mathgroup@wolfram.com
  • Delivered-to: mathgroup-newout@smc.vnet.net
  • Delivered-to: mathgroup-newsend@smc.vnet.net

Bellow is a program designed to track a[t] Rm[t] and T[t]. I want to use it to find out by what amount a[t]changes over the course of the cycle shown. I want to track from the initial point where Rm[t] = 1 to the final point when it again returns to 1. Ideally I would like to plot this change in a[t] as a function of R to view how the asymmetry changes with R. From there I would really like to plot a 3D array of a[t] as a function of R as well as alpha. Thanks ahead of time for any help. Bellow the initial syntax is a program I have started to try and accomplish this.

M = 100;
g = 1;
\[Alpha] = 1;
ti = 0;(*initial time*)tf = 12000;(*final plot time*)m = 1;
RB[t] = (m*T[t]^3*g)/(2*\[Pi])^(3/2)*(m/T[t])^(3/2)*E^(-m/T[t]);
R = .0002;
V = 5*10^(-12)*M^4;
i = Sqrt[(g*m^4 - V + g/(2*\[Pi])^(3/2)*E^-1)/(3*M^2)]*m;
rR = g*T[t]^4;(*Energy Density-Radiation*)pR = 
 g*T[t]^4;(*Pressure-Radiation*)pM = 0;(*Pressure-Phi Particle*)r = 
 rR + Rm[t] - V;(*Overall Energy Density*)p = 
 pR + V;(*Overall Pressure*)s1 = 
 NDSolve[{a''[t] == -a[t]*1/(6*M^2) (2*g*T[t]^4 + 2 V + Rm[t]), 
   T'[t] == (R*Rm[t])/(4*
        T[t]^3) + (Sqrt[3/2]*\[Alpha]^2)/(4*Sqrt[T[t]]*
         m^(11/2))*((RB[t])^2 - Rm[t]^2) - (a'[t]*T[t])/a[t], 
   Rm'[t] == -R*
      Rm[t] + (Sqrt[3/2]*\[Alpha]^2*T[t]^(5/2))/
       m^(11/2)*((RB[t])^2 - Rm[t]^2) - 3 a'[t]/a[t]*Rm[t], T[0] == m,
    Rm[0] == g/(2*\[Pi])^(3/2)*E^-1, a'[0] == i, a[0] == 1}, {a[t], 
   T[t], Rm[t]}, {t, ti, tf}, Method -> "BDF"]

Plot[Evaluate[Rm[t]/RB[t] /. s1], {t, ti, tf}, PlotRange -> Automatic,
  AxesLabel -> {Time, Rm[t]}]
Plot[Evaluate[T[t] /. s1], {t, ti, tf}, PlotRange -> Automatic, 
 AxesLabel -> {Time, T[t]}]
Plot[Evaluate[a[t] /. s1], {t, ti, tf}, PlotRange -> Automatic, 
 AxesLabel -> {Time, a[t]}]
ParametricPlot[Evaluate[{T[t], a[t]} /. s1], {t, ti, tf}, 
 PlotRange -> Automatic]



curveColor[n_Integer] := 
  Hue[Mod[2/3 - 2 + Sqrt[5]*(n - 1), 1], 0.6, 0.6];

Manipulate[
 Module[{M, g, ti, tf, m, RB, V, i, rR, pR, pM, r, p, imgSize, 
   aspRatio, s1, plt}, M = 100;
  g = 1;
  ti = 0;(*initial time*)tf = 12000;(*final plot time*)m = 1;
  RB[t] = (m*T[t]^3*g)/(2*\[Pi])^(3/2)*(m/T[t])^(3/2)*E^(-m/T[t]);
  V = 5*10^-12*M^4;
  i = Sqrt[(g*m^4 - V + g/(2*\[Pi])^(3/2)*E^-1)/(3*M^2)]*m;
  rR = g*T[t]^4;(*Energy Density-Radiation*)pR = g*T[t]^4;(*Pressure-
  Radiation*)pM = 0;(*Pressure-Phi Particle*)
  r = rR + Rm[t] - V;(*Overall Energy Density*)
  p = pR + V;(*Overall Pressure*)imgSize = 150;
  s1 = NDSolve[{a''[t] == -a[t]*1/(6*M^2) (2*g*T[t]^4 + 2 V + Rm[t]), 
      T'[t] == (R*Rm[t])/(4*
           T[t]^3) + (Sqrt[3/2]*\[Alpha]^2)/(4*Sqrt[T[t]]*
            m^(11/2))*((RB[t])^2 - Rm[t]^2) - (a'[t]*T[t])/a[t], 
      Rm'[t] == -R*
         Rm[t] + (Sqrt[3/2]*\[Alpha]^2*T[t]^(5/2))/
          m^(11/2)*((RB[t])^2 - Rm[t]^2) - 3 a'[t]/a[t]*Rm[t], 
      T[0] == m, Rm[0] == g/(2*\[Pi])^(3/2)*E^-1, a'[0] == i, 
      a[0] == 1}, {a[t], T[t], Rm[t]}, {t, ti, tf}, 
     Method -> "BDF"][[1]];
  Row[{Show[
     plt = Plot[
       Evaluate[{Tooltip[Rm[t]/RB[t], "Rm[t]"], 
          Tooltip[20 T[t], "T[t]"], Tooltip[a[t], "a[t]"]} /. s1], {t,
         ti, tf}, Frame -> True, Axes -> False, ImageSize -> imgSize, 
       PlotRange -> {-1, 22}, AspectRatio -> 2.5], 
     FrameTicks -> {{Automatic, (FrameTicks /. 
            AbsoluteOptions[plt, FrameTicks])[[2]] /. {y_, 
           yl_?NumericQ, r__} :> {y, yl/20, 
           r}}, {(FrameTicks /. 
            AbsoluteOptions[plt, FrameTicks])[[1]] /. {x_, 
           xl_?NumericQ, r__} :> {x, xl/1000, r}, Automatic}}, 
     FrameLabel -> {{Row[{Style["Rm[t]", Bold, curveColor[1]], 
          "   &   ", Style["a[t]", Bold, curveColor[3]]}], 
        Style["T[t]", Bold, curveColor[2]]}, {"Time/1000", None}}], 
    ParametricPlot[Evaluate[{T[t], a[t]} /. s1], {t, ti, tf}, 
     Frame -> True, FrameLabel -> {"T[t]", "a[t]"}, Axes -> False, 
     AspectRatio -> 2, ImageSize -> imgSize, 
     PlotRange -> {{0.14, 0.5}, {2, 7}}, 
     ColorFunction -> (ColorData["Rainbow"][#3] &)]}]], {{\[Alpha], 0,
    "|\[Alpha]|"}, 0, 5, Appearance -> "Labeled"}, {{R, 0}, 0, 0.001, 
  Appearance -> "Labeled"}]



Here is my attempt.

M = 100;
g = 1;
\[Alpha] = 1;
ti = 0;(*initial time*)tf = 12000;(*final plot time*)m = 1;
RB[t] = (m*T[t]^3*g)/(2*\[Pi])^(3/2)*(m/T[t])^(3/2)*E^(-m/T[t]);
R = .0002;
V = 5*10^(-12)*M^4;
i = Sqrt[(g*m^4 - V + g/(2*\[Pi])^(3/2)*E^-1)/(3*M^2)]*m;
rR = g*T[t]^4;(*Energy Density-Radiation*)pR = 
 g*T[t]^4;(*Pressure-Radiation*)pM = 0;(*Pressure-Phi Particle*)r = 
 rR + Rm[t] - V;(*Overall Energy Density*)p = 
 pR + V;(*Overall Pressure*)s1 = 
 NDSolve[{a''[t] == -a[t]*1/(6*M^2) (2*g*T[t]^4 + 2 V + Rm[t]), 
   T'[t] == (R*Rm[t])/(4*
        T[t]^3) + (Sqrt[3/2]*\[Alpha]^2)/(4*Sqrt[T[t]]*
         m^(11/2))*((RB[t])^2 - Rm[t]^2) - (a'[t]*T[t])/a[t], 
   Rm'[t] == -R*
      Rm[t] + (Sqrt[3/2]*\[Alpha]^2*T[t]^(5/2))/
       m^(11/2)*((RB[t])^2 - Rm[t]^2) - 3 a'[t]/a[t]*Rm[t], T[0] == m,
    Rm[0] == g/(2*\[Pi])^(3/2)*E^-1, a'[0] == i, a[0] == 1}, {a[t], 
   T[t], Rm[t]}, {t, ti, tf}, Method -> "BDF"]

Plot[Evaluate[Rm[t]/RB[t] /. s1], {t, ti, tf}, PlotRange -> Automatic,
  AxesLabel -> {Time, Rm[t]}]
Plot[Evaluate[T[t] /. s1], {t, ti, tf}, PlotRange -> Automatic, 
 AxesLabel -> {Time, T[t]}]
Plot[Evaluate[a[t] /. s1], {t, ti, tf}, PlotRange -> Automatic, 
 AxesLabel -> {Time, a[t]}]
ParametricPlot[Evaluate[{T[t], a[t]} /. s1], {t, ti, tf}, 
 PlotRange -> Automatic]



  • Prev by Date: Re: Crashing every other launch?
  • Next by Date: Re: Crashing every other launch?
  • Previous by thread: Tracking Differential Equations
  • Next by thread: maximization with array of constraints