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]