MathGroup Archive 2012

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

Search the Archive

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]
>>



  • 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