# Problems with doing calculus on Interpolating functions

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!

• 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
– Rojo
Feb 22, 2013 at 18:17
• Hiya, sorry about that. I've condensed the problem down a bit...thanks for pointing that out! Feb 22, 2013 at 18:28
• You need to provide numerical values for all the constants too. Feb 22, 2013 at 18:54
• done! sorry about that too... Feb 22, 2013 at 18:55

## 1 Answer

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


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

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