# How can I apply calculus to functions obtained from NDSolve?

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?

• No one can usefully comment if you don't show any code. – Ymareth Aug 21 '14 at 12:48
• [Gamma] = 6; g = -9.8; NDSolve[{y''[t] + [Gamma]*(y'[t])^2 == g, q''[t] == -[Gamma]*(q'[t])^2, y[0] == 0, q[0] == 0, y'[0] == 1.5, q'[0] == 7}, {q, y}, {t, 0, 10}] – Jason Aug 21 '14 at 12:57
• Sorry i don't know how to convert them to looks better here's the answers: {{q-> InterpolatinFunction[{{0.,0.31767}},<>],y-> InterpolatinFunction[{{0.,0.31767}},<>]}} – Jason Aug 21 '14 at 12:59
• You're using NDSolve which gives a numerical answer, essentially a table of values of the q or gamma co-ordinate against time wrapped up as an InterpolatingFunction. If you want the functional form of the line you either need to use DSolve (which may not be able to solve this) or fit a function to the result from NDSolve. – Ymareth Aug 21 '14 at 13:01
• This question has received three reasonable answers based on different interpretations of what is being asked. The OP has not indicated that any of these interpretations are helpful to him; therefore, I have voted to close this question as "not clear what you are asking". – m_goldberg Aug 22 '14 at 3:34

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

• Michael E2 , sorry i didn't get it first, but i get it now, thank you very much for your answer – Jason Aug 22 '14 at 10:42
• In routine maintenance I was preparing to delete this old closed question. If you want this answer to be preserved I suggest editing this abandoned question to better fit this answer. – Mr.Wizard Jul 25 '15 at 13:56

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.

• There is also FindFit that typically can do the job within Mma without attractind additional programs. – Alexei Boulbitch Aug 21 '14 at 14:52
• Thanks Alexei Boulbitch just a question of you , how can i extract my points? because it's giving me this as i said :q-> InterpolatinFunction[{{0.,0.31767}},<>],y-> InterpolatinFunction[{{0.,0.31767}},<>]}} as an answer how can i extract points from this? – Jason Aug 21 '14 at 15:02
• I agree that FindFit will work, but the problem with it is that the user has to specify a functional form. There may be little reason here to think that the function is linear or some simple polynomial. The virtue of the approaches I outlined is that the user merely specifies primitives to be used in the evolution of a functional form. It's kind of a difference between a pure "statistics" approach and more of a "machine learning" approach. Both have virtues. – Seth Chandler Aug 21 '14 at 15:36

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"}]]


• thank you for your answer but i'm not looking for path graph i already have it!, i'm looking for the approximate function cz i'm gonna work on that function, i need some information that i can't get them from the graph! – Jason Aug 21 '14 at 16:46
• This isn't about the path graph. I only show it to provide visual confirmation that f` is the expected path function and, therefore, the approximate function you are seeking. – m_goldberg Aug 21 '14 at 23:30
• I have edited my answer to clarify the point I want to make. – m_goldberg Aug 22 '14 at 0:09
• I got it Thank you – Jason Aug 22 '14 at 8:56