# Numerical Integration of 2nd order ODE

I need to integrate an equation of the form

$f''(x) g_1+f'(x) g_2=a g_3$

where $a$ is a constant and $g_1,g_2,g_3$ are functions of $x$ whose values on the nodes of integration are known. The boundary conditions are $f(x=0)=0$, $f(x=1)=f_0$, $f_0$ is a constant. $g_1,g_2,g_3$ are arrays of data, a value for each node.

My question is: should I interpolate $g_1,g_2,g_3$ and then use the NDSolve command, or it is possible to employ directly a numerical tecnique?

• Is this question related to the Software Wolfram Mathematica or more of a question on mathematical / numerical solution methods? If so, then you may want to try math.stackexchange.com or scicomp.stackexchange.com Nov 29 '16 at 13:54
• I believe that this problem can be solved using Mathematica, so it is relevant with this topic.
– DK13
Nov 29 '16 at 14:00
• How is mathematica supposed to differentiate between a g3 and b g4? There's not enough information here. Nov 29 '16 at 14:03
• I think the answer is yes, you need to interpolate first if you want to use NDSolve. If you want real help, please provide a minimal working example of the constants, gs and the boundary conditions to solve the problem. Nov 29 '16 at 14:05
• @DK13 I wouldn't bother trying to re-invent the NDSolve Nov 29 '16 at 16:18

a finite difference soluton:

dx = 1/8;
grid = Range[-1., 1., dx];
g1 = Table[1. + x^2, {x, grid}];
g2 = Table[x, {x, grid}];
g3 = Table[Cos[x], {x, grid}];
a = 2;
unk = Array[f, {Length@grid}];

result = Transpose[{grid,
unk /. Solve[Join[{unk[] == 0, unk[[-1]] == -2.5},
Table[
(f[i + 1] - 2 f[i] + f[i - 1])/dx^2 g1[[i]] +
(f[i + 1] - f[i-1])/2/dx g2[[i]] ==  a g3[[i]],
{i, 2, Length[grid] - 1}]], unk][]}];
ListPlot[result, Joined -> True, PlotMarkers -> Automatic] here is a check: fpp g1 + fp g2 and a g3 plotted together.

fp = Derivative[Interpolation[result]]
fpp = Derivative[Interpolation[result]]
Show[{
ListPlot[
Table[ {grid[[i]],
fpp[grid[[i]]] g1[[i]] + fp[grid[[i]]] g2[[i]] }, {i,
Length@grid}], Joined -> True],
ListPlot[Table[ {grid[[i]], a g3[[i]] }, {i, Length@grid}]]}] note this formulation does not enforce the equation at the ends.

• Thank you! it was really helpful...
– DK13
Nov 29 '16 at 20:55

Here's a finite element method way:

(* data for OP's problem *)
grid = Range[-1., 1., 1./8];
g1 = Table[1. + x^2, {x, grid}];
g2 = Table[x, {x, grid}];
g3 = Table[Cos[x], {x, grid}];

(* solution to OP's problem *)
Needs["NDSolveFEM"];
emesh = ToElementMesh[
"Coordinates" -> List /@ grid,
"MeshElements" -> {LineElement[Partition[Range@Length@grid, 2, 1]]}];
G1 = Interpolation[Transpose[{grid, g1}], InterpolationOrder -> 1];
G2 = Interpolation[Transpose[{grid, g2}], InterpolationOrder -> 1];
G3 = Interpolation[Transpose[{grid, g3}], InterpolationOrder -> 1];
Block[{a = 2},
{sol} =
NDSolve[{f''[x] G1[x] + f'[x] G2[x] == a G3[x],
DirichletCondition[f[x] == 0, x == -1]},
f, {x} ∈ emesh]
];

(* for comparison *)
Block[{a = 2},
{sol2} =
NDSolve[{f''[x] (1 + x^2) + f'[x] x == a Cos[x], f[-1] == 0,
f'[-1] == (f'[-1] /. sol)}, f, {x, -1, 1}]
];

Plot[{f[x] /. sol, f[x] /. sol2}, {x, -1, 1}] • it is really great, but where is the Dirichlet condition at $x=1$ ??
– DK13
Nov 29 '16 at 20:11
• @DK13 When I first read your question it didn't have any BCs....I left your original Q on my browser and didn't refresh and see the update until just now. Anyway it's simple enough to add it: NDSolve[{..., DirichletCondition[f[x] == 0, x == 0], DirichletCondition[f[x] == f0, x == 1]},...], isn't it? Nov 29 '16 at 21:35
• I did it already! Thank you very much!!!! @Michael E2
– DK13
Nov 30 '16 at 11:42