# How do I numerically integrate a numerical solution of a differential equation? [duplicate]

Suppose I have received a numerical solution to a differential equation via

ODE = f''[x] + f[x] - 1/f[x]^2 == 0;
sol = NDSolve[{ODE, f == 1, f' == 2}, f[x], {x, 0, 10}]


Now, I would like to use the solution further

1) for numerical integration

g[t] := NInt[f[t] /. sol, x, {x, 0, t}]


So that, e.g.

g


would give me a value according to the above definition of g[t]. The above does not work.

2) also for more complicated evaluations, such as

h[t] := Exp[NInt[f[t] /. sol, x, {x, 0, t}]]


whereby I'd like to plot the data and store into a file afterwards, i.e. something like:

Plot[h[t],{t, 0, 10}]


which does not work either.

I know these are basic operations but so far I was unable to get through the obstacles of coding the above.

• This might be helpful: reference.wolfram.com/language/howto/… – yohbs May 10 '17 at 17:16
• For 1, you can still use NDSolve[]: NDSolveValue[{f''[x] + f[x] - 1/f[x]^2 == 0, f == 1, f' == 2, g'[x] == f[x], g == 0}, {f, g}, {x, 0, 10}]. – J. M.'s ennui May 10 '17 at 22:38

There seem to be some mistakes in your approach. First I start with your ODE and then I plot the solution to check that it is sound.

ODE = f''[x] + f[x] - 1/f[x]^2 == 0;
sol = NDSolve[{ODE, f == 1, f' == 2}, f[x], {x, 0, 10}];
Plot[Evaluate[f[x] /. First@sol], {x, 0, 10}] Note the Evaluate this insures that the replacement is done once before plotting begins.

Your next line of code is

h[t] := Exp[NInt[f[t] /. sol, x, {x, 0, t}]]


What is NInt? This does not seem to be a Mathematica function. If you mean NIntegrate then this does not have the form you wrote. Do you mean

g[t_] := NIntegrate[Evaluate[f[x] /. First@sol], {x, 0, t}]


Note that I have made a pattern of g[t_] and removed the first x in the expression. Is this what you wanted?

A better approach would be to include this integration in the original ODE. So starting again we have

ODE1 = {f''[x] + f[x] - 1/f[x]^2 == 0, g'[x] == f[x], f == 1,
f' == 2, g == 0};
sol = NDSolve[ODE1, {f[x], g[x]}, {x, 0, 10}];


The function for g is now available

Plot[Evaluate[g[x] /. sol], {x, 0, 10}] There are several syntax mistakes. The correct way to define the function g[t] is the following:

g[t_] := NIntegrate[f[x] /. sol[], {x, 0, t}]

g
(* 8.65985 *)

1. In order to define a function with a variable you need to use a pattern t_ which means that every appearance of t at the other side of the assignation operator will be substituted by the input value to the function. Writing g[t] assigns an expression to this particular combination of symbols, but does not create a function.

2. Notice the correct syntax for NIntegrate. You can consult it in Mathematica executing ?NIntegrate or using the help manuals.

3. I have written /.sol[] instead of /.sol because NDSolve (and many other functions in Mathematica) returns a list of replacement rules, although in this case is a list with only one item. If you use sol you will get a list of all the results of the replacement (again, in this case a list of only one element, but a list). Of course, there's no problem with that and may be your desired result. In other case, you should avoid extracting the desired solution from the list.

Method I

ODE = f''[x] + f[x] - 1/f[x]^2 == 0;
{f[x_]} = {f[x]} /.First@NDSolve[{ODE, f == 1, f' == 2}, f[x], {x, 0, 10}];
g[t_] := NIntegrate[f[x], {x, 0, t}];
g
h[t_] := Exp[g[t]];
Plot[h[t], {t, 0, 10}, PlotRange -> All]


Method II

ODE = f''[x] + f[x] - 1/f[x]^2 == 0;
{fsol} = NDSolveValue[{ODE, f == 1, f' == 2}, {f}, {x, 0, 10}];
g[t_] := NIntegrate[fsol[x], {x, 0, t}]
h[t_] := Exp[g[t]]
Plot[h[t], {t, 0, 10}, PlotRange -> All] 