Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: How to relate two functions

• To: mathgroup at smc.vnet.net
• Subject: [mg127476] Re: How to relate two functions
• From: Bob Hanlon <hanlonr357 at gmail.com>
• Date: Sat, 28 Jul 2012 02:38:39 -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
• References: <20120727085736.04021684A@smc.vnet.net>

M = 100;
g = 1;
\[Alpha] = 1;
ti = 0;(*initial time*)
tf = 700;(*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 = 5*10^-4;
V = 5*10^(-12)*M^4;
i = Sqrt[(g*m^4 - V)/(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])/T[t]^3 -
       (Sqrt[3/2]*\[Alpha]^2)/(Sqrt[T[t]]*m^(11/2))*
        (((g*T[t]^3*m)/(2*\[Pi])^(3/2)*
             (m/T[t])^(3/2)*
             E^(-m/T[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)*
        ((g/(2*\[Pi])^(3/2)*(m/T[t])^(3/2)*
             E^(-m/T[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]];

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

Find values of t for which T[t] == 1

r = FindRoot[T[t] == 1 /. s1,
      {t, #}] & /@ {575, 1} // Chop // Quiet

{{t -> 633.688}, {t -> 0}}

Corresponding values of a[t]

val = a[t] /. s1 /. r

{1.74697, 1.}

Difference

Subtract @@ %

0.74697

ParametricPlot[
 {T[t], a[t]} /. s1, {t, ti, tf},
 PlotRange -> Automatic,
 Frame -> True,
 FrameLabel -> {"T[t]", "a[t]"},
 Axes -> False,
 Epilog -> {Red,
   AbsoluteDashing[{5, 5}],
   Line[Prepend[{1, #} & /@ val, {0.6, val[[1]]}]]}]


Bob Hanlon


On Fri, Jul 27, 2012 at 4:57 AM, William Duhe <williamduhe at hotmail.com> wro=
te:
> So here I have several functions plotted, a[t], B[t],Rm[t]. What I need i=
s syntax that will allow me to see the value of a[t] when T[t] is 1. This s=
hould be easy to see that at first when T[t] is 1 a[t] is also one, but as =
the function turns around T[t] becomes 1 again, at about t = 575, and a[t=
] has risen by a certain amount. I need to quantify this amount.I have plot=
ted a parametric plot of a[t] vs T[t] and what I need to do is find the Y-I=
NTERCEPTS of that plot in order to get the data I need, but other more soph=
isticated methods would also be welcome. Thanks ahead of time for any help!
>
> The syntax is provided bellow.
>
>
> M = 100;
> g = 1;
> \[Alpha] = 1;
> ti = 0;(*initial time*)
> tf = 700; (*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 = 0.0005;
> V = 5*10^(-12)*M^4;
> i = Sqrt[(g*m^4 - V)/(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])/
>      T[t]^3 - (Sqrt[3/2]*\[Alpha]^2)/(
>       Sqrt[T[t]]*m^(
>        11/2))*(((g*T[t]^3*m)/(2*\[Pi])^(3/2)*(m/T[t])^(3/2)*
>           E^(-m/T[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)*((g/(2*\[Pi])^(3/2)*(m/T[t])^(3/2)*E^(-m/T[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]}]
> Plot[{Evaluate[a[t] /. s1], Evaluate[T[t] /. s1]}, {t, ti, tf},
>  PlotRange -> Automatic, AxesLabel -> {Time, a[t]}]
> ParametricPlot[Evaluate[{T[t], a[t]} /. s1], {t, ti, tf},
>  PlotRange -> Automatic]
>

• References:
  • How to relate two functions
    • From: William Duhe <williamduhe@hotmail.com>