Mathematica Asked by haozz on February 1, 2021
This is a question about the fluid mechanics equation, which is solved by a similarity solution ($f(t)$, here).
I’m trying to solve the following boundary value problem with shooting method (taken from $(2)(3)(4)$ of this paper):
$f(t)-t f^{prime}(t)+aleft(f(t)^{3} f^{prime prime prime}(t)right)^{prime}=0$
$f(0)=1, f^{prime}(0)=f^{prime prime prime}(0)=0, f^{prime prime}(infty)=0, f^{prime}(infty)=1$
Five boundary conditions are given, in order to determine the unknown parameter $a$.
I choose the ParametricNDSolveValue
with the first four boundary conditions, the fifth condition is used conducting the shooting method. Infinity is replaced by t==100000
, but there are some errors with the results:
pfun = ParametricNDSolveValue[{f[t] == t f'[t] - a D[f[t]^3 f'''[t], t],
f[0] == 1, f'[0] == f'''[0] == 0, f''[100000] == 0},
f'[100000], {t, 0, 100000}, {a}]
FindRoot[pfun[a] - 1, {a, 2}]
Unfortunately, Mathematica gives something like these:
Power::infy: Infinite expression 1/0.^3 encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.
To sum up, my questions are: how can I figure out this ODE to check whether the boundary condition at infinity (in my shooting algorithm I take infinity as t = 100 000
) is satisfied? Is my setting wrong?
Thanks!
Update:
When I set the xi as x=10, it still doesn’t work. There is a singular at t=0. Errors are shown as:
Power::infy: Infinite expression 1/0.^3 encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
ParametricNDSolveValue::ndnum: Encountered non-numerical value for a derivative at t$3391 == 0.`.
However, when I change ‘a’ to ‘-a’, it seems to get a strange answer, which is beyond my expectation. In fact, the value of ‘a’ should be around 1.22, as stated in an article.
Update2:
The final purpose is to fix this equation:
$f-x f^{prime}+a left(f^{R+2}left|f^{prime prime prime}right|^{R-1} f^{prime prime prime}right)^{prime}=0$
$f(0)=1, f^{prime}(0)=f^{prime prime prime}(0)=0, f^{prime prime}(infty)=0, f^{prime}(infty)=1$
Find ‘a’ for a specific value of ‘R’, the prior question is under the condition R=1. I have tried as:
R = 2;
{fsol, asol} =
NDSolveValue[{f[t] ==
t f'[t] -
a[t] D[f[t]^(R + 2) (Abs [f'''[t]])^(R - 1)*f'''[t], t],
a'[t] == 0, f[0] == 1, f'[0] == f'''[0] == 0, f''[10] == 0,
f'[10] == 1}, {f, a}, {t, 0, 10}];
Plot[{fsol[t], asol[t]}, {t, 0, 10}]
y1 = asol[1]
if R=1, y1 = 1.3417, which is corresponding to answer of @xzczd;
When R takes other values, errors appear:
Power::infy: Infinite expression 1/0. encountered.
NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.`.
So this problem may be difficult to solve,owing to the singular at t==0.
Setting Shooting method outside of *NDSolve*
with FindRoot
helps, this might be related to the arguable backslide of "Shooting"
method.
inf = 10;
{eq, bc} = {f[t] == t f'[t] - a D[f[t]^3 f'''[t], t], {f[0] == 1, f'[0] == 0,
f'''[0] == 0, f''[inf] == 0, f'[inf] == 1}};
pfun = ParametricNDSolveValue[{eq, bc[[;; 3]], f''[0] == c}, f, {t, 0, inf}, {a, c}]
parasol = FindRoot[{pfun[a, c]'[inf] == 1, pfun[a, c]''[inf] == 0}, {{a, 2}, {c, 2}},
MaxIterations -> 500]
(* {a -> 1.3417, c -> 0.632144} *)
The result is slightly different from the a = 0.818809^-1
mentioned in the paper, but it's actually a better one, at least when $infty$ is approximated as 10
:
Block[{a = 0.818809^-1},
pfuntst = ParametricNDSolveValue[{eq, bc[[;; 3]], f''[0] == c}, f, {t, 0, inf}, {c}]]
parasoltst = FindRoot[{pfuntst[c]'[inf] == 1}, {{c, 2}}]
(* {c -> 0.661846} *)
{pfuntst[c]'[inf], pfuntst[c]''[inf]} /. parasoltst
(* {1., -0.0041385} *)
{pfun[a, c]'[inf], pfun[a, c]''[inf]} /. parasol
(* {1., 4.93118*10^-15} *)
As you can see, my $f''(infty)$ is closer to $0$.
The 2 solutions are quite close to each other by the way:
{pfuntst[c] /. parasoltst, pfun[a, c] /. parasol} // ListLinePlot
You may adjust the parameters to e.g. inf = 5
to check the result further.
The code above is tested on v11.3, v12.0.1 and v12.1.1. In v9.0.1 FindRoot
spits out a njnum
warning and returns unevaluated, which seems to be a bug.
Answered by xzczd on February 1, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP