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;

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;
> pR =
> 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;
>> 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]
>>

```

• Prev by Date: Re: How to relate two functions
• Next by Date: Re: Obtaining densest x% of a Bivariate Distribution
• Previous by thread: Re: How to relate two functions
• Next by thread: Re: How to relate two functions