How can I apply calculus to functions obtained from NDSolve?

Question

Originally, I asked the question below, but the real underlying issue is as follows: When we solve an ODE numerically, I get the answer like this:

{{y-> InterpolatingFunction[{{0.,0.386145}},<>]}}

How can I use this result to do some calculus on the approximate function (for example, differentiate, optimize, integrate, or even plot)?

Here is example code that produces such functions:

γ = 6; 
g = -9.8; 
NDSolve[{y''[t] + γ*(y'[t])^2 == g, q''[t] == -γ*(q'[t])^2, 
  y[0] == 0, q[0] == 0, 
  y'[0] == 1.5, q'[0] == 7}, {q, y}, {t, 0, 10}]

---

Original question: I solved two ODEs, which are a function of t, numerically. The first ODE is the vertical equation of motion and the second one is the horizontal equation of that motion. Then I tried to find the path equation. I used ParametricPlot to plot the path. Now I want the equation of that path. How can I express it using the solutions of my ODEs?

Answer 1

From the comments, we can set up the OP's DEs as follows and show they can be solved exactly.

First the system is the direct product of two independent systems, so let's separate them.

γ = 6; g = -98/10;
yIVP = {y''[t] + γ*(y'[t])^2 == g, y[0] == 0, y'[0] == 15/10};
qIVP = {q''[t] == -γ*(q'[t])^2, q[0] == 0, q'[0] == 7};

nysol0 = NDSolve[yIVP, {y}, {t, 0, 10}]
nqsol0 = NDSolve[qIVP, {q}, {t, 0, 10}]

NDSolve::ndsz: At t == 0.3176700784294118`, step size is effectively zero; singularity or stiff system suspected. >>

(*
  {{y -> InterpolatingFunction[{{0., 0.31767}}, <>]}}
  {{q -> InterpolatingFunction[{{0., 10.}}, <>]}}
*)

We can try to solve the systems exactly, and DSolve quickly returns:

dysol0 = DSolve[yIVP, y, t]
dqsol0 = DSolve[qIVP, q, t]

Solve::incnst: Inconsistent or redundant transcendental equation. After reduction, the bad equation is -196+466 Cos[7 Sqrt[30] C[1]]^2 == 0. >>

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

(* dysol0:
  {{y -> Function[{t}, 
      1/6 (-I π - Log[7 Sqrt[2/233]] + 
         Log[Cos[7 Sqrt[6/5] (t + 1/7 Sqrt[5/6] ArcCos[-7 Sqrt[2/233]])]])]},
   {y -> Function[{t}, 
    1/6 (-Log[7 Sqrt[2/233]] + 
       Log[Cos[7 Sqrt[6/5] (t - 1/7 Sqrt[5/6] ArcCos[7 Sqrt[2/233]])]])]}}
*)

DSolve::bvnul: For some branches of the general solution, the given boundary conditions lead to an empty solution. >> (* dqsol0: {} *)

OK, so we've got some work to do. Let's try for general solutions (again, DSolve returns quickly):

dysol = DSolve[First[yIVP], y, t]
dqsol = DSolve[First[qIVP], q, t]
(*
  {{y -> Function[{t}, C[2] + 1/6 Log[Cos[7 Sqrt[6/5] (t - 5 C[1])]]]}}
  {{q -> Function[{t}, C[2] + 1/6 Log[6 t - C[1]]]}}
*)

That's encouraging. Let's investigate further.

Note that Rest[yIVP] /. First[dysol] gives the initial conditions:

Rest[yIVP] /. First[dysol]
Rest[qIVP] /. First[dqsol]
(*
{C[2] + 1/6 Log[Cos[7 Sqrt[30] C[1]]] == 0, (7 Tan[7 Sqrt[30] C[1]])/Sqrt[30] == 3/2}

{C[2] + 1/6 Log[-C[1]] == 0, -(1/C[1]) == 7}
*)

They don't look that bad, but if we try Solve, we get a result similar to DSolve:

Solve[Rest[yIVP] /. First[dysol], {C[1], C[2]}]
Solve[Rest[qIVP] /. First[dqsol], {C[1], C[2]}]

Solve::incnst: Inconsistent or redundant transcendental equation. After reduction, the bad equation is -196+466 Cos[7 Sqrt[30] C[1]]^2 == 0. >>

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

{{C[1] -> -(ArcCos[-7 Sqrt[2/233]]/(7 Sqrt[30])), 
  C[2] -> 1/6 (-I π - Log[7 Sqrt[2/233]])},
 {C[1] -> ArcCos[7 Sqrt[2/233]]/(7 Sqrt[30]), 
  C[2] -> -(1/6) Log[7 Sqrt[2/233]]}}

{}

Let's try Reduce instead:

Reduce[Rest[yIVP] /. First[dysol], {C[1], C[2]}]
Reduce[Rest[qIVP] /. First[dqsol], {C[1], C[2]}]
(*
  C[3] ∈ Integers && 
   C[1] == (ArcTan[(3 Sqrt[15/2])/7] + π C[3])/(7 Sqrt[30]) && 
   C[2] == -(1/6) Log[Cos[7 Sqrt[30] C[1]]]

  C[1] == -(1/7) && C[2] == Log[7]/6
*)

Ah, looks like success! The solution for y needs a choice for C[3] but it doesn't matter what integer we pick because its effect on C[1], which appears inside Cos, makes no difference in the answer. Here then are the solutions:

ycoeff = Reduce[Rest[yIVP] /. First[dysol], {C[1], C[2]}] /. C[3] -> 0 // ToRules
qcoeff = Reduce[Rest[qIVP] /. First[dqsol], {C[1], C[2]}] // ToRules
(*
  {C[1] -> (π + ArcTan[(3 Sqrt[15/2])/7])/(7 Sqrt[30]), 
   C[2] -> -(1/6) Log[Cos[7 Sqrt[30] C[1]]]}

  {C[1] -> -(1/7), C[2] -> Log[7]/6}
*)

dysol0 = dysol //. ycoeff
dqsol0 = dqsol /. qcoeff
(*
  {{y -> Function[{t}, -(1/6) Log[Cos[7 Sqrt[30] C[1]]] + 
      1/6 Log[Cos[7 Sqrt[6/5] (t - (5 (π + ArcTan[(3 Sqrt[15/2])/7]))/(7 Sqrt[30]))]]]}}

  {{q -> Function[{t}, Log[7]/6 + 1/6 Log[6 t - -(1/7)]]}}
*)

As a check, we'll plot our symbolic solutions on top of the numeric solutions:

t1 = 0;
t2 = t /. First@Solve[
    Cos[7 Sqrt[6/5] (t - 1/7 Sqrt[5/6] ArcCos[7 Sqrt[2/233]])] == 0 &&0 < t < 1/2, t];

Plot[{y[t] /. First@nysol0, y[t] /. First@dysol0}, {t, t1, t2}]
Plot[{q[t] /. First@nqsol0, q[t] /. First@dqsol0}, {t, 0, 10}]

Looks good. The OP mentioned DSolve being slow. That is true on the original system. Solving each component separately is much faster.

Update

Re: How can i Approximate the function to do some calculus on it?

The numerical (interpolating) functions nysol0, nqsol0 returned by NDSolve is an approximation of the exact functions dysol0, dqsol0 that you can do calculus on. You can also do calculus on the exact functions.

Examples:

We scale q to better match the range of its derivative q'.

Plot[{10 q[t], q'[t]} /. First@nqsol0 // Evaluate, {t, t1, t2}]

If using First@<> is inconvenient, extract the InterpolatingFunction.

y0 = y /. First@nysol0; (* y0 is now a function *)
{{t1, t2}} = y0["Domain"];

FindRoot[y0'[t] == 0, {t, t2/2}]
(* {t -> 0.112822} *)

FindMaximum[y0[t], {t, t2/2, t1, t2}]
(* {0.0721726, {t -> 0.112822}} *)

Integrate[y0[t], {t, t1, t2}]
(* -0.00328963 *)

Answer 2

There are at least two approaches by which you could obtain a closed form emulation of your answer. Both involve extracting a list of points that are part of the solution. 1. Inexpensive but takes you outside Mathematica. Export the list of points to a CSV file. Obtain (free trial) a program called Eureqa (http://www.nutonian.com/products/eureqa/) that uses genetic programming to evolve functions that will closely emulate the points you have selected. It's pretty easy to use this product.
2. Expensive. Get a package called DataModeler from EvolvedAnalytics (http://www.evolved-analytics.com/) that runs within Mathematica and likewise uses genetric programming to evolve functions that will closely emulate the points you have selected. This package is extremely powerful and advanced. Might be overkill if this is some sort of homework or single hobby problem.

Perhaps others will have better solutions. There is, of course, LinearModelFit on some nth degree polynomial, but there you are basically pre-selecting the functional form and not letting the data speak for itself.

Answer 3

You can get a path function directly from the solution of your ODEs.

I don't understand your ODE's, so I'm going to work with a much simpler system, which gives the path of particle moving under constant gravity in a vacuum.

g = -9.8;
numSoln = 
  NDSolve[{
       y''[t] == g, q''[t] == 0., 
       y[0] == 0., q[0] == 0, 
       y'[0] == 50., q'[0] == 15.}, 
    {q, y}, {t, 0., 10.}];

From numSoln, the numerical solution of the ODEs, the approximate path function is given by

f[t_] = {q[t], y[t]} /. numSoln;

The above expression for f is independent of the precise form of the ODEs, so it is applicable to your problem as well. The following plot is made to visually confirms the result.

ParametricPlot[f[t], {t, 0., 10.}, AxesLabel -> {q, y}]

The components q and y as functions of t can also be recovered from f.

Plot[Evaluate @ f[t], {t, 0., 10.},
  AxesLabel -> {t, None},
  PlotLegends -> SwatchLegend[Automatic, {"q", "y"}]]