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$
4
  • 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, 2013 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, 2013 at 18:28
  • 1
    $\begingroup$ You need to provide numerical values for all the constants too. $\endgroup$
    – b.gates.you.know.what
    Feb 22, 2013 at 18:54
  • $\begingroup$ done! sorry about that too... $\endgroup$
    – Gokotai
    Feb 22, 2013 at 18:55

1 Answer 1

Reset to default
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$
3
  • $\begingroup$ as a side note, you could also use NDSolveValue that would avoid that ... /. ScaleFactor[[1]] business. $\endgroup$
    – user21
    Feb 23, 2013 at 8:45
  • $\begingroup$ @ruebenko Not me :) ... still on v.8 $\endgroup$
    – Dr. belisarius
    Feb 23, 2013 at 8:48
  • $\begingroup$ ah, OK, then as an idea for the OP, and maybe a small incentive to upgrade ;-) $\endgroup$
    – user21
    Feb 23, 2013 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.