4
$\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$
4
$\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

|improve this answer|||||
$\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

|improve this answer|||||
$\endgroup$

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.