Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

First off, appologies for what may sound like a newbie question, as I am very new to using Mathematica.

I am trying to find a way to get Mathematica to give me an expression that would describe the data I get from my code, but I'm not sure how. I've played around with FindFit and can't seem to get it to work. Also, I have an integral whose integrand relies on two variables, but I'd only like to integrate with respect to one and keep the other constant. I can do this at the "defining a function" stage but am unsure as to how to display the information thereafter. Any help would be much appreciated! Here is my code

a := 3.24077*10^-20 (* km \[Rule] Mpc *)
b := 3.16888*10^-14 (* s \[Rule] MYear *)
c := a/b*(2.99792*10^5) (*Mpc/ MYear*)
H0 := a/b*71 (*1/MYear*)
G := a^3/b^2*6.67398*10^-20 (*Mpc^3/(Kg*MYear^2)*)
\[CapitalOmega]M := 0.27
\[CapitalOmega]\[Gamma] := 8.24*10^-4
\[CapitalOmega]\[CapitalLambda] := 0.73
f := 1*10^-4
y := 100
\[Lambda] := (1/(a^2)) 1*(10^-58)(*1/Mpc^2*)
H[z_] := H0 Sqrt[\[CapitalOmega]M (1 + 
      z)^3 + \[CapitalOmega]\[Gamma] (1 + 
      z)^4 + \[CapitalOmega]\[CapitalLambda]]
\[Rho]crit := 3/(8 \[Pi]*G)*(H0)^2 
\[Rho]Pert[r_] := \[CapitalOmega]M*\[Rho]crit*(1 + f*Exp[-(r/y)^2])
M1[r_] := 4 \[Pi] Integrate[\[Rho]Pert[r]*r^2, r]
E1[r_] := 
 1/2*((H[999]*r)/(c*(1000)))^2 - ((G*M1[r]*(1000))/(
   c^2*(r))) - (\[Lambda]*r^2)/(3 (1000)^2)
ScaleFactor = 
  NDSolveValue[
   SetPrecision[{(Sqrt[R[r, t]] (D[R[r, t], t])) == Sqrt[
      2*G*M1[r] + 2*c^2*E1[r]*R[r, t] + (2 c^2 \[Lambda]*R[r, t]^3)/
       3], R[r, 0] == r/1000}, 200], R, {r, 1, 20000}, {t, 1, 15000}];
Plot3D[Re[ScaleFactor[r, t]], {r, 1, 20000}, {t, 1, 15000}, 
 PlotRange -> Automatic, Axes -> True, 
 AxesLabel -> {"Comoving Radius (Mpc)", "Time (MegaYears)"}, 
 LabelStyle -> Directive[Bold, Black]]
FindFit[ScaleFactor[r, 
  t], \[Alpha]*r^n + \[Beta]*t^m, {\[Alpha], \[Beta], n, m}, {r, t}]
ParticleHorizon[r_, t_] := 
 c*Integrate[Sqrt[1 + 2*E1[r]]/D[ScaleFactor[r, t], r], t]
share|improve this question
    
As currently posted, this strikes me as a "can you debug my code" question. –  m_goldberg Mar 11 '13 at 14:20

1 Answer 1

up vote 0 down vote accepted

Maybe

sample=Join @@ Table[{r, t, ScaleFactor[r, t]}, {r, 1, 20000, 200}, {t, 1, 15000,  150}];
FindFit[sample, \[Alpha]*r^n + \[Beta]*t^m, {\[Alpha], \[Beta], n, m}, {r, t}]
(* {\[Alpha] -> 0.00118915, \[Beta] -> 0.00206976, n -> 1.93354, m -> 1.93373} *)

?

share|improve this answer
    
Thanks alot, that works great for the FindFit! However, I still get an error when I try to do a ParallelTable for ParticleHorizon, which is defined as follows: ParticleHorizondata = ParallelTable[ ParticleHorizon[r, t], {r, 1, 10}, {t, 1, 10}] –  Gokotai Mar 11 '13 at 14:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.