-1
$\begingroup$

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?

$\endgroup$
11
  • $\begingroup$ 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 $\endgroup$ Nov 29 '16 at 13:54
  • $\begingroup$ I believe that this problem can be solved using Mathematica, so it is relevant with this topic. $\endgroup$
    – DK13
    Nov 29 '16 at 14:00
  • $\begingroup$ How is mathematica supposed to differentiate between a g3 and b g4? There's not enough information here. $\endgroup$
    – Feyre
    Nov 29 '16 at 14:03
  • $\begingroup$ 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. $\endgroup$ Nov 29 '16 at 14:05
  • 3
    $\begingroup$ @DK13 I wouldn't bother trying to re-invent the NDSolve $\endgroup$
    – Chris K
    Nov 29 '16 at 16:18
4
$\begingroup$

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[[1]] == 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][[1]]}];
ListPlot[result, Joined -> True, PlotMarkers -> Automatic]

enter image description here

here is a check: fpp g1 + fp g2 and a g3 plotted together.

fp = Derivative[1][Interpolation[result]]
fpp = Derivative[2][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}]]}]

enter image description here

note this formulation does not enforce the equation at the ends.

$\endgroup$
1
  • $\begingroup$ Thank you! it was really helpful... $\endgroup$
    – DK13
    Nov 29 '16 at 20:55
4
$\begingroup$

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["NDSolve`FEM`"];
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}]

Mathematica graphics

$\endgroup$
3
  • $\begingroup$ it is really great, but where is the Dirichlet condition at $x=1$ ?? $\endgroup$
    – DK13
    Nov 29 '16 at 20:11
  • $\begingroup$ @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? $\endgroup$
    – Michael E2
    Nov 29 '16 at 21:35
  • $\begingroup$ I did it already! Thank you very much!!!! @Michael E2 $\endgroup$
    – DK13
    Nov 30 '16 at 11:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.