Re: How to relate two functions
- To: mathgroup at smc.vnet.net
- Subject: [mg127481] Re: How to relate two functions
- From: Bob Hanlon <hanlonr357 at gmail.com>
- Date: Sat, 28 Jul 2012 02:40:20 -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>
Manipulate[
Module[
{M, g, \[Alpha], ti, tf, m, RB, V, i, rR,
pR, pM, r, p, s1, t1, val},
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]);
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]];
t1 = FindRoot[T[t] == 1 /. s1, {t, 600}];
Column[{
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 -> {{-5, 705}, {0, 3.75}},
AxesLabel -> {"Time", "functions"},
ImageSize -> 350,
Exclusions -> Automatic],
t1 = Quiet[
Chop[FindRoot[T[t] == 1 /. s1, {t, #}] & /@ {600, 1}]],
val = a[t] /. s1 /. t1,
Subtract @@ val}]],
{{R, 5.*^-4}, 1.*^-6, 0.007, Appearance -> "Labeled"}]
Bob Hanlon
On Fri, Jul 27, 2012 at 10:32 AM, Bob Hanlon <hanlonr357 at gmail.com> wrote:
> 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> w=
rote:
>> So here I have several functions plotted, a[t], B[t],Rm[t]. What I need =
is syntax that will allow me to see the value of a[t] when T[t] is 1. This =
should 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 plo=
tted a parametric plot of a[t] vs T[t] and what I need to do is find the Y-=
INTERCEPTS of that plot in order to get the data I need, but other more sop=
histicated 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>
- How to relate two functions