1
$\begingroup$

The following 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)*)

\[Rho]crit1 := 3/(8 \[Pi]*G)*(H0)^2 
    M1[r_] := (4 \[Pi] )/3*\[Rho]crit1*(r)^3 
    E1[r_] := 0
    ScaleFactor = 
      NDSolve[{Sqrt[R[r, t]] (D[R[r, t], t]) == 
          Sqrt[(2*G*M1[r])/c + 2 c*E1[r]*R[r, t]], R[r, 0] == r/1000}, 
        R, {r, 1, 10000}, {t, 1, 10000}] // FullSimplify;
    Plot3D[Evaluate[R[r, t] /. %], {r, 1, 10000}, {t, 1, 10000}]
    ParticleHorizon[r_, g_] := 
     NIntegrate[(c Sqrt[1 + 2 E1[r]])/
      D[Evaluate[R[r1, t] /. ScaleFactor, r1 -> r], r], {t, 0, g}]
    Plot3D[ParticleHorizon[r,g], {r, 1, 10}, {g, 1, 10}]

My problem is that once I've solved the differential equation to get R[r,t] (which it solves without a problem, giving me an interpolating function), I have trouble differentiating the result with respect to r, to be used in my ParticleHorizon function. I think it might be because the differential operator uses DSolve or NDSolve intrinsically.

Here, a,b,c,f,G and lambda are constants and M1[r] is a boundary condition that we need to solve the differential equation for R[r,t]. E1[r] has been set to 0 for now for simplicity.

Any ideas on how to proceed?

Also, the rule to change r1->r was just me trying to see if it works, which it doesn't...

Thanks in advance!

|improve this question
$\endgroup$
  • 1
    $\begingroup$ Hey, welcome! Try to find a minimal example that illustrate your problem, so that it shows your effort, is less work for the people trying to help you, and is more general and helps future visitors more $\endgroup$ – Rojo Feb 22 '13 at 18:17
  • $\begingroup$ Hiya, sorry about that. I've condensed the problem down a bit...thanks for pointing that out! $\endgroup$ – Gokotai Feb 22 '13 at 18:28
  • 1
    $\begingroup$ You need to provide numerical values for all the constants too. $\endgroup$ – b.gates.you.know.what Feb 22 '13 at 18:54
  • $\begingroup$ done! sorry about that too... $\endgroup$ – Gokotai Feb 22 '13 at 18:55
3
$\begingroup$

You could do:

ParticleHorizon[r_, g_] := 
   NIntegrate[(c Sqrt[1 + 2 E1[r]])/Evaluate[D[R[r, t] /. ScaleFactor[[1]]], r], {t, 0, g}]

ListPlot3D[
        Flatten[Table[{r, g, ParticleHorizon[r, g]}, {r, 1, 10}, {g, 1, 10}], 1], 
        ColorFunction -> "SouthwestColors", Mesh -> True, 
        MeshFunctions -> {#3 &}, 
        PlotStyle -> Directive[Opacity[0.8], Specularity[White, 50]]]

Mathematica graphics

I used ListPlot3D just because it is faster than Plot3D[]

|improve this answer
$\endgroup$
  • $\begingroup$ as a side note, you could also use NDSolveValue that would avoid that ... /. ScaleFactor[[1]] business. $\endgroup$ – user21 Feb 23 '13 at 8:45
  • $\begingroup$ @ruebenko Not me :) ... still on v.8 $\endgroup$ – Dr. belisarius Feb 23 '13 at 8:48
  • $\begingroup$ ah, OK, then as an idea for the OP, and maybe a small incentive to upgrade ;-) $\endgroup$ – user21 Feb 23 '13 at 8:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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