# Interpolating a ParametricFunction

Original domain of my independent variable t is 0 to 10 nanometers, but in Mathematica (as in most numerical software), it's good idea to solve differential equation in the more 'usual' unit scale 0 to 10.

pfun = ParametricNDSolveValue[ {x''[
t] + a x'[t] + b^2 x[t] == 0, x[0] == 1,
x'[0] == 0}, x, {t, 0, 10}, {a,b}]


How can I scale back t in the resulting ParametricFunction to nanometers?

Note: With NDSolve we used before version 9 (which introduced ParametricNDSolveValue), the result was an InterpolatingFunction, which was easily scaled back by using

Interpolation[{#*10^(-9), pfun[#]} & /@Range[0, 10, 10^(-1)]]

-
There seems to me be an extra factor of 10 in NDSolve example. If not, my answer is off by a factor of 10. –  Michael E2 Aug 19 '13 at 17:49
The expression {#*10*10^(-9), pfun[#]} was allright, since 10*10^(-9) equals to 10 nm or 1.*^-9 in more compact notation. Therefore I edited your answer to make it compatible. –  Cendo Aug 20 '13 at 8:09
Are you sure? 10*10^(-9) == 1.*^-8 returns True on my machine. –  Michael E2 Aug 20 '13 at 11:23
Or to put it another way, if Range goes from 0 to 10, then the expression #*10*10^(-9) goes from 0 to 10*10*10^(-9) or 100 * 10^-9. I assumed you were converting between nm & m (but maybe that's wrong?). I thought maybe the 10 in Range or the factor of 10 might have been meant to be a 1 or left out. –  Michael E2 Aug 20 '13 at 13:24
Sorry, my stupid mistake with #*10*10^(-9). Now I corrected it in the answer. And yes, i was converting nm &m. –  Cendo Aug 20 '13 at 13:40

In the differential equation x is the dependent variable, a function of t and which function depends on parameters a, b. So we have x given by pfun[a, b][t], and I take it from the example, that t is to be rescaled. If so, then this seems serviceable (for a = -2, b = 1):

xfn = pfun[-2, 1][Rescale[#, {0, 10^-9}]] &


Example:

xfn[4.*^-9]
(* -163.794 *)


If the output value x is to be rescaled then

xfn = Rescale[pfun[-2, 1][#], {0, 10^-9}] &


and

xfn[4.]
(* -1.63794 * 10^11 *)


Generally, one might define

x[a_?NumericQ, b_?NumericQ, t_?NumericQ] := pfun[a, b][Rescale[t, {0, 10^-9}]]


Update

Alternatively, note that pfun[a, b] yields an InterpolatingFunction, which you can use as you are used to (although I realize that might be inconvenient in many applications). There is a drawback to the way the function is re-interpolated in the OP that is worth addressing.

The InterpolatingFunction has more data, and data that is more precise, at certain points in the domain than at the equally spaced sample points given by Range[0, 10, 10^(-1)]]. NDSolve stores the values of the function and its first and second derivatives at input values called a "Grid".

Here's a way to rescale the InterpolatingFunction returned by ParametricNDSolveValue. We have to scale the derivatives by different factors. (You can justify it either by the "Chain Rule" or by dimensional analysis of $dx/dt$ and $d^2x/dt^2$.)

If ifn is an InterpolatingFunction1, you can get the grid of input numbers and the corresponding function values with

ifn["Grid"]
ifn["ValuesOnGrid"]

(* {{x1}, {x2}, ... }
{ y1,   y2,  ... } *)


and ifn[{"Grid", "ValuesOnGrid"}] returns both in a list; ifn["Coordinates"] also returns the input numbers in a flat list. To get the values of the derivative (AFAIK), you need to evaluate ifn' and ifn'' on the grid ("Coordinates" is convenient here).

rescale = {1*^-9, 1, 1*^9, 1*^18};
solfun[a_?NumericQ, b_?NumericQ] := solfun[a, b] =
Interpolation[Transpose[
rescale * Join[
pfun[a, b][{"Grid", "ValuesOnGrid"}],
Through[{pfun[a, b]', pfun[a, b]''}[First[pfun[a, b]["Coordinates"]]]]
]
]]


We can compare it to a re-interpolation on a regular grid.

ifun = Interpolation[{# * 10^(-9), pfun[-2, 1][#]} & /@ Range[0, 10, 10^(-1)]];


The mean relative error of solfun is much better:

Table[Abs[(ifun[t] - pfun[-2, 1][Rescale[t, {0, 10^-9}]])/
pfun[-2, 1][Rescale[t, {0, 10^-9}]]], {t, 0, 10^-8, 10^-13}] // Mean
Table[Abs[(solfun[-2, 1][t] - pfun[-2, 1][Rescale[t, {0, 10^-9}]])/
pfun[-2, 1][Rescale[t, {0, 10^-9}]]], {t, 0, 10^-8, 10^-13}] // Mean

(* 9.03955*10^-6
6.50966*10^-15 *)


The maximum relative error is likewise several orders of magnitude better for solfun. The horizontal axis is just about at the level of Log10[\$MachineEpsilon]. (There's something funny going on at 10^-9.)

ListPlot[Log10@
{Table[Abs[(ifun[t] - pfun[-2, 1][Rescale[t, {0, 10^-9}]])/
pfun[-2, 1][Rescale[t, {0, 10^-9}]]],
{t, 0, 10^-8, 10^-13}],
Table[Abs[(solfun[-2, 1][t] - pfun[-2, 1][Rescale[t, {0, 10^-9}]])/
pfun[-2, 1][Rescale[t, {0, 10^-9}]]],
{t, 0, 10^-8, 10^-13}]},
DataRange -> {0, 10^-8}]


Ref.: See this answer for more about InterpolatingFunctions.

-
Actually, I searched for a solution which would also scale t inside pfun before knowing the values of a,b (mainly for purpose of speeding up subsequent solutions of ParametricNDSolveValue. I don't know if that is possible to do, since it certainly has to do with how Mathematica calculates solution internally. –  Cendo Aug 20 '13 at 8:36
As far as I can tell, ParametricNDSolveValue after checking of arguments, sets up a family of DEs. When, say, pfun[-2, 1] is first called, the DE is solved, via NDSolve probably, and stores ("caches") the InterpolatingFunction solution, which is reused on subsequent calls. So for instance Plot[pfun[-2, 1][t], {t, 0, 10}] (DE solved once) is about 16 times faster per function call than Plot3D[pfun[a, b][10.], {a, 0, 1}, {b, 0, 1}] (DE solved many times, once for each a, b). Btw, Rescale might be fast enough. –  Michael E2 Aug 20 '13 at 11:46
For me your Plot[] example is 60 (sixty) times faster than Plot3D (measured with Timing[]). –  Cendo Aug 20 '13 at 13:44
I had about 100 times faster in all, but Plot had 580 function calls (to pfun) and Plot3D had 3928 function calls. For brevity, I didn't include the option EvaluationMonitor :> (cnt++) in my previous comment. Set cnt = 0 before the plot. –  Michael E2 Aug 20 '13 at 13:59

If I understand your question right, you only need to put correct initial conditions. Indeed, your equation

eq = x''[t] + a x'[t] + b^2 x[t] == 0


is linear. I understand that you mean that x is measured in meters and t in seconds. If I am right, let us make a substitution: x[t] -> 10^m*y[t] with m=9

eq /. x -> (10^m*y[#] &)
Map[Expand@Divide[#, 10^m] &, %]


The result is:

(*
10^m b^2 y[t] + 10^m a Derivative[1][y][t] +
10^m (y^\[Prime]\[Prime])[t] == 0

b^2 y[t] + a Derivative[1][y][t] + (y^\[Prime]\[Prime])[t] == 0
*)


So we get the same equation for any m, but now for y=y[t]. This is the general result due to the linearity of the equation in question. Now, if m=9, y is measured in nanometers. You only need to set the initial conditions also in nanometers.

By the way, did you really mean that x[0]==1 meter? Probably you had in mind 1 nanometer? Anyway, it seems more convenient to look for the analytical solution:

 DSolve[{y''[t] + a y'[t] + b^2 y[t] == 0, y[0] == 1, y'[0] == 0}, y, t]

(*   {{y -> Function[{t}, (1/(
2 Sqrt[a^2 - 4 b^2]))(-a E^(1/2 (-a - Sqrt[a^2 - 4 b^2]) t) +
Sqrt[a^2 - 4 b^2] E^(1/2 (-a - Sqrt[a^2 - 4 b^2]) t) +
a E^(1/2 (-a + Sqrt[a^2 - 4 b^2]) t) +
Sqrt[a^2 - 4 b^2] E^(1/2 (-a + Sqrt[a^2 - 4 b^2]) t))]}}
*)

-
I just corrected my question that I meant both t and x in nanometers. Please edit your answer appropriately –  Cendo Aug 20 '13 at 8:45
@Cendo Your question is still obscure. I still do not see, where the problem lies. All your unknowns and parameters should have correct dimensions that agree to one another. That is dimensions of both a and b are 1/nm and x[0]==1 nm. You need not to do anything else. Just solve it in Mathematica, and write down dimensions to the right of the solution. Done. –  Alexei Boulbitch Aug 21 '13 at 7:43