# 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...

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

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 '13 at 8:45
@ruebenko Not me :) ... still on v.8 – Dr. belisarius Feb 23 '13 at 8:48
ah, OK, then as an idea for the OP, and maybe a small incentive to upgrade ;-) – user21 Feb 23 '13 at 8:52