Interpolating a ParametricFunction

Question

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

Answer 1

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.

Answer 2

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))]}}
    *)