5
$\begingroup$

Question 101668 was not particularly clear and was put on hold for that reason. The author did clarify it a bit, but it seems unlikely that the question will be reopened. Through an exchange of comments with the author, I believe that I understand what he asking, which I present here and for which I offer an answer. Perhaps others can offer better answers.

The OP wishes to solve the system of equations

{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0, g''[x] + p[[i]] f[x] g'[x] == 0,
    f[0] == f'[0] == 0, g'[0] == -1, f'[5] == 2, g[5] == 0}

for a set of parameter values.

p = {0.01, 0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.4, 0.6, 0.72, 1, 5, 7, 10, 100, 1000};

The obvious approach is

a = ParametricNDSolve[{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0, 
    g''[x] + p0 f[x] g'[x] == 0, f[0] == f'[0] == 0, g'[0] == -1, 
    f'[5] == 2, g[5] == 0}, {f, g}, {x, 0, 5}, {p0}];

over the domain {x, 0, 5}, which works well for p as large as 7

Plot[Evaluate[{g[7][x] /. a, f[7][x] /. a}], {x, 0, 5}, PlotRange -> All]

enter image description here

but appears to become unstable for larger values of p. (According to the OP, Maple can solve this problem for large p without difficulty.) One might hope that adding the options

Method -> "StiffnessSwitching", MaxSteps -> Infinity

would provide large p answers, but it is exceedingly slow. Setting

WorkingPrecision-> 30

instead helps a bit, but not much. It too is slow. The question, then, is - how to obtain solutions for large p in a reasonable amount of time.

$\endgroup$
0
5
$\begingroup$

One approach to obtaining large p results involves two observations. First, f is independent of p and can be solved at the outset.

sf = First[f /. NDSolve[{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0, 
    f[0] == f'[0] == 0, f'[5] == 2}, f, {x, 0, 5}]];

Second, a first integral exists for the ODE involving g.

lgp = -Integrate[sf[x], x];

where lgp is Log[g'[x]] for p = 1. Thus, two applications of Integrate produce a solution.

hp[p0_?NumericQ] := FunctionInterpolation[-Exp[p0 lgp], {x, 0, 5}, 
    InterpolationPoints -> 1000, InterpolationOrder -> 2]
h[p0_?NumericQ, x_?NumericQ] := -Integrate[hp[p0][x0], {x0, x, 5}]

Note that InterpolationPoints must be fairly large to provide reasonable accuracy at large p.

Plot[Evaluate[h[1000, x]], {x, 0, 5}, PlotRange -> All, AxesLabel -> {x, g[1000]}]

enter image description here

$\endgroup$
0
$\begingroup$

This problem should be related to this one. Based on the experience got in that post, first analyse the equation a bit:

DSolve[{g''[x] + p0 f[x] g'[x] == 0(*, g[0] == 0*), g'[0] == 0}, g[x], x]
{{g[x] -> C[2]}}

Apparently the solution will be {{g[x] -> 0}} if g[0] == 0 (though DSolve will fail if this b.c. is directly included in it), which suggests that we can force g[x] and g'[x] to be exact 0 once they're close enough to 0 to avoid the error accumulation. So here comes the solution:

sf = NDSolveValue[{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0, f[0] == f'[0] == 0, 
    f'[5] == 2}, f, {x, 0, 5}];

sg = With[{p0 = 1000, f = sf, e = $MachineEpsilon}, 
  NDSolveValue[{g''[x] + p0 f[x] g'[x] == 0, g'[0] == -1, g[5] == 0, 
    WhenEvent[{g[x] < e, g'[x] < e}, {g[x] -> 0, g'[x] -> 0}]}, g, {x, 0, 5}]]

Plot[sg[x], {x, 0, 5}]

enter image description here

$\endgroup$

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .