Trouble with the shooting method for boundary value problem of a 4th-order ODE

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)+a\left(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. • Hello, welcome to Mathematica.SE. Then, it should be f[t]== rather than f[t]=. Also, you may want to read this post to learn how to format the code properly. Sep 8, 2020 at 5:25 • This looks like a nonlinear eigenvalue problem to me, with$a$as the eigenvalue. It is fourth order and$f(t) \equiv 0 $is a solution, so you should only need to specify 4 boundary conditions, not 5. I don't know much at all about nonlinear eigenvalue problems, but for a linear eigenvalue problem, the 5th condition would let you scale the eigenfunction once you have found an eigenvalue. Generally setting infinity to be 100000 is much too large for these kind of problems, try something like 10 initially. Sep 8, 2020 at 5:29 • Getting down to your error message, your equation is singular at$t=0\$, and you'll have to perturb it off that point. Sep 8, 2020 at 5:42
• Thank you for the answers, but it has always failed. Sep 8, 2020 at 7:00
• The updated statements have been added. Sep 8, 2020 at 7:51

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.

• Sir, thanks for your answers. The value of 'a' mentioned in the paper is about 1/0.818809 = 1.22, the final result 1.34 seems to be great. But, as I added the new equation above, when an absolute value and the changed exponent 'R' are added, Can this equation still be solved? We want to find the relationship between the exponent 'R' and the undecided 'a'. Sep 10, 2020 at 4:24
• @haozz What have you tried? Sep 10, 2020 at 5:20
• when R=1, we have used another method to get the value of 'a', as a=1.3417. When R takes other values, errors appear. See the updated contents above. Sep 10, 2020 at 6:15