Using an interpolating solution of NDSolve into an algebraic expression

Question

I would like to use the interpolating solution of NDSolve into an algebraic expression without creating a system of a differential equation and algebraic equation.

The example is the following, the differential eq. is

NDSolve[{
  x'[t] == -Sqrt[(1/(4 - 2 Sqrt[2]) + (4 - 2 Sqrt[2])/(2 x[t]))^2 - 1],
  x[0] == 4}, 
  x, {t, 0, 10}]

The algebraic expression is

m[t] Sqrt[1 + x'[t]^2] - m[t]^2/(2 x[t]) + 1 == 0

that I need to solve for m[t]

Answer 1

You have your NDSolve

xsoln = NDSolve[{x'[
     t] == -Sqrt[(1/(4 - 2 Sqrt[2]) + (4 - 2 Sqrt[2])/(2 x[t]))^2 - 
       1], x[0] == 4}, x, {t, 0, 10}]

which returns a replacement rule. This replacement rule will work with derivatives, just try Plot[{x[t], x'[t]} /. xsoln, {t, 0, 10}, Evaluated -> True] to see this.

Then you solve the algebraic equation and act upon it with this replacement rule, generating a new set of replacement rules

msoln = Solve[m[t] Sqrt[1 + x'[t]^2] - m[t]^2/(2 x[t]) + 1 == 0, 
   m[t]] /. First@xsoln

You can plot the two solutions using this replacement rule,

Plot[m[t] /. msoln, {t, 0, 10}, Evaluated -> True]

Answer 2

Initialization

One alternative is to construct new interpolating functions from the solutions.

dsol = x /. 
   First@NDSolve[{x'[
        t] == -Sqrt[(1/(4 - 2 Sqrt[2]) + (4 - 
                2 Sqrt[2])/(2 x[t]))^2 - 1], x[0] == 4}, 
     x, {t, 0, 10}];

msol = Solve[m[t] Sqrt[1 + x'[t]^2] - m[t]^2/(2 x[t]) + 1 == 0, m[t]]
(*
  {{m[t] -> x[t] Sqrt[1 + Derivative[1][x][t]^2] - 
      Sqrt[x[t]] Sqrt[2 + x[t] + x[t] Derivative[1][x][t]^2]},
   {m[t] ->  x[t] Sqrt[1 + Derivative[1][x][t]^2] + 
      Sqrt[x[t]] Sqrt[2 + x[t] + x[t] Derivative[1][x][t]^2]}}
*)

Solution

See below for discussion.

removearg = h_[t] :> h;
vals = {
   x -> dsol["ValuesOnGrid"],
   x' -> dsol'["ValuesOnGrid"],
   x'' -> NDSolve`FiniteDifferenceDerivative[Derivative[2], 
     Flatten@dsol["Grid"], dsol["ValuesOnGrid"], "DifferenceOrder" -> 2]};

mfns = MapThread[
   Interpolation@Transpose@{dsol["Grid"], ##} &,
   {m[t] /. msol /. removearg /. vals,
    D[m[t] /. msol, t] /. removearg /. vals}
   ]

Plot[Through[mfns[t]], {t, 0, 10}, Evaluated -> True]

Discussion

To get an order-3 Hermite interpolation, like the one produced by NDSolve, we need to supply the function and derivative values of m on an interpolation grid. Here is a fake example just to let you inspect the pieces of the computation without being overwhelmed by the volume of data in dsol.

MapThread[
 Transpose@{List /@ Range[3], ##} &,
 {m[t] /. msol /. h_[t] :> h /. {
    x -> Range[3]^3,
    x' -> 3. Range[3]^2},
  D[m[t] /. msol, t] /. h_[t] :> h /. {
    x -> Range[3]^3,
    x' -> 3. Range[3]^2,
    x'' -> 6. Range[3]}}
 ]
(*
  {{{{1}, -0.301824, 0.456501}, {{2}, -0.0830097, 0.0823127}, {{3}, -0.0370107, 0.0246378}},
   {{{1}, 6.62638, 29.9014}, {{2}, 192.749, 480.253}, {{3}, 1459.04, 2430.31}}}
*)

Since m[t] contains x'[t], then m'[t] will contain x''[t]. There are a couple of problems with x''[t]. In cubic Hermite interpolation, the second derivative is not guaranteed to be continuous. The discontinuities will be on the interpolation grid, so estimating x''[t] at the grid points involves a choice. The second problem is that the initial steps taken by NDSolve are usually quite small. The purpose seems to be to take small steps while enough steps are accumulated so that the finite difference derivatives are reasonably accurate and the error estimate can be trusted. One can see that the second derivative oscillates a bit wildly at the beginning. It over a very small interval, so its importance depends on the importance of the behavior of m[t] for t near 0. (See Fig. 1.)

Fig. 1. A plot of x''[t] showing the oscillation and discontinuities near t == 0. The range of the oscillations is in fact larger than 0.3, or about 1000 times as great as shown.

Plot[dsol''[t], {t, 0, 1/2}, PlotStyle -> Thick]

We can use NDSolve`FiniteDifferenceDerivative to estimate x''[t] as seen in the code for vals. Since NDSolve computes and stores x'[t] in the solution, these values on the grid points may be assumed to be accurate. So we might estimate x''[t]`` from eitherdsol[t]ordsol'[t]. Each way is roughly equivalent and better than usingdsol''[t]. I chose usingdsol[t]` above.

Fig. 2. Estimates of x''[t]: dsol''[t] (blue), the (finite difference) first derivative of dsol'[t] (yellow), and the (finite difference) second derivative of dsol[t] (green).

ListLinePlot[
 {dsol''["ValuesOnGrid"],
  NDSolve`FiniteDifferenceDerivative[Derivative[1], 
   Flatten@dsol["Grid"], dsol'["ValuesOnGrid"], 
   "DifferenceOrder" -> 2],
  NDSolve`FiniteDifferenceDerivative[Derivative[2], 
   Flatten@dsol["Grid"], dsol["ValuesOnGrid"], 
   "DifferenceOrder" -> 2]},
 PlotRange -> All]

As a final note, the rule removearg = h_[t] :> h used in the code converts m[t] to m, x[t] to x, x'[t] to x' and so forth.