# How to solve this fourth-order singular nonlinear ODE?

Hi I would like to hear your suggestions on solving the following fourth-order singular nonlinear ODE regarding $u=u(x)$

$\alpha u'''' + u'u''' + (u'')^2 = 0$

where prime denotes derivative w.r.t. x and $\alpha$ is a small quantity (therefore singular). The above ODE is subject to the boundary conditions

$u(1)=1, u(-1)=0, u''(1)=-10, u'''(-1)=0$.

The tricky part in the ODE is the boundary conditions. They are asymmetric and peculiar. Do you have any suggestions? Thanks a lot.

• Note that $u'u'''+u''^2=(u'u'')'=\frac12(u'^2)''$ which allows to reduce the order of the ODE. – Lutz Lehmann Jun 13 '16 at 17:19
• Hi Lutz, I know this transformation. And I indeed also reached the general solution of the ODE, but the problem is the boundary conditions. I don't know how to use them to determine the integration constants appearing in the general solution, see below for Mariusz's procedure. – jengmge Jun 13 '16 at 18:46

$\alpha u''''(x)+u''(x)^2+u'''(x) u'(x)=0$

Integrating once:

$\alpha u'''(x)+u'(x) u''(x)=\text{c1}$

Integrating again:

$\alpha u''(x)+\frac{1}{2} u'(x)^2=\text{c1} x+\text{c2}$

substitution: $u'(x)=v(x)$,now you have differential equations of the first order:

$\alpha \frac{\mathrm{d} v(x)}{\mathrm{d} x}+\frac{v(x)^2}{2}=\text{c1} x+\text{c2}$

solving:

$v(x)=\frac{2^{2/3} \alpha \sqrt[3]{\frac{\text{c1}}{\alpha ^2}} \left(\text{c3} \text{Ai}'\left(\frac{(\text{c1} x-\text{C2}) \sqrt[3]{\frac{\text{c1}}{\alpha ^2}}}{\sqrt[3]{2} \text{c1}}\right)+\text{Bi}'\left(\frac{(\text{c1} x-\text{C2}) \sqrt[3]{\frac{\text{c1}}{\alpha ^2}}}{\sqrt[3]{2} \text{c1}}\right)\right)}{\text{c3} \text{Ai}\left(\frac{(\text{c1} x-\text{C2}) \sqrt[3]{\frac{\text{c1}}{\alpha ^2}}}{\sqrt[3]{2} \text{c1}}\right)+\text{Bi}\left(\frac{(\text{c1} x-\text{C2}) \sqrt[3]{\frac{\text{c1}}{\alpha ^2}}}{\sqrt[3]{2} \text{c1}}\right)}$

back of substitution and integrating:

$u(x)=2 \alpha \ln \left(\text{c3} \text{Ai}\left(\frac{\text{c1} x-\text{c2}}{\sqrt[3]{2} \left(\frac{\text{c1}}{\alpha }\right)^{2/3} \alpha }\right)+\text{Bi}\left(\frac{\text{c1} x-\text{c2}}{\sqrt[3]{2} \left(\frac{\text{c1}}{\alpha }\right)^{2/3} \alpha }\right)\right)+\text{c4}$

Where $\text{Ai}$ and $\text{Bi}$ is Airy function.

To find the integrations constants c1,c2,c3,c4 You need to solve 4 transcedental equations.

Edited:

The only way to find integrations constants its only numerically.

ClearAll["Global*"]
Remove["Global*"]

sol = First@DSolve[\[Alpha]*u''''[x] + u'[x] u'''[x] + u''[x]^2 == 0, u[x], x];
U = u[x] /. sol /. C[1] -> c1 /. C[2] -> c2 /. C[3] -> c3 /. C[4] -> c4;
eq = {1 == U /. x -> 1, 0 == U /. x -> -1, -10 == D[U, {x, 2}] /. x -> 1, 0 == D[U, {x, 3}] /. x -> -1};

\[Alpha] = 1;

sol2 = FindRoot[eq, {{c1, 1}, {c2, 2}, {c3, 1}, {c4, 1}}]


$\{\text{c1}\to -42.3064-29.7117 i,\text{c2}\to -30.2477-39.8327 i,\text{c3}\to 0.000396293\, -0.998705 i,\text{c4}\to 0.786101\, +1.89902 i\}$

Check the boundary conditions:

 Re[U /. sol2] /. x -> 1,                 (* u[1] == 1 *)
Re[U /. sol2] /. x -> -1,                (* u[-1] == 0 *)
Re[D[(U /. sol2), {x, 2}]] /. x -> 1,    (* u''[1] == -10 *)
Re[D[(U /. sol2), {x, 3}]] /. x -> -1    (*u'''[-1] == 0 *)

(* {1., -1.5099*10^-13, -10., -1.31024*10^-10} *)

Plot[{Re[U /. sol2], Evaluate[Re[D[(U /. sol2), {x, 2}]]],
Evaluate[Re[D[(U /. sol2), {x, 3}]]]}, {x, -1, 1},
PlotLegends -> {"u[x]", "u''[x]", "u'''[x]"},
PlotRange -> {Automatic, {-11, 2}}]


Edited:

for $\alpha = 1/50$

\[Alpha] = 1/50;
sol2 = FindRoot[eq, {{c1, -10.7}, {c2, -11.88}, {c3, 0}, {c4, 1.09}},
Method -> "Newton", MaxIterations -> 1000000]


Checking the boundary conditions:

 Re[U /. sol2] /. x -> 1,                 (* u[1] == 1 *)
Re[U /. sol2] /. x -> -1,                (* u[-1] == 0 *)
Re[D[(U /. sol2), {x, 2}]] /. x -> 1,    (* u''[1] == -10 *)
Re[D[(U /. sol2), {x, 3}]] /. x -> -1    (*u'''[-1] == 0 *)

(* {0.999742, 0.64927, -10.102, 0.27074} *)


The results are not accurate.The tricky parts is a find good starting points in FindRoot.

• Hi Mariusz, indeed, I also arrived at similar results as you, but I expressed the results in terms of the Bessel functions. But it is the task of determining these integration constants that drives me crazy, that's the trickiest part as the boundary conditions are so peculiar. I used Mathematica to do that, but it was no possible to solve them. Do you have any suggestions or tricks to do that? Thanks. – jengmge Jun 13 '16 at 18:42
• Thanks a lot. I did try Mathematica with the Bessel function and the integration constants, then the things are going crazy. It's a good idea just to solve the equation directly. I will check later. Thanks. – jengmge Jun 14 '16 at 7:46