Problem solving Third order non-linear differential equation in Mathematica

Question

I am trying to find an analytical solution of the following 3rd order non-linear differential equation in Mathematica: $a (f'(x))^2+f'''(x)=0$ with boundary conditions $f(0)=0$, $f'(0)=0$, $f(1)=1$, $f''(1)=1$

This how I input it in Mathematica.

DSolve[{a (f'[x])^2+f'''[x]==0, f[0]=0, f'[0]==0,f[1]=1,f''[1]=0},f[x],x]

... Dsolve: For some branches of the general solution, unable to solve for the conditions

    {}

Would you give any suggestions how to make it work?

Answer 1

Elliptic functions are not optimally implemented in Mathematica, there are many examples on this site affirming this statement, and it would be reasonable to emphasise what does not work perfectly. This is the reason why DSolve does not provide an exact solution with coefficients determined. However we can choose a bit more pedestriaan approach.

Taking a look at the boundary conditions we can see that the coefficient $a$ cannot be equal to $0$. Let's denote $f'(x)=g(x)$. Now we have $g''+a\;g^2=0$. We can multiply the last equation by $g'$ and integrate it once to get ${g'}^2-4(-\frac{a}{6})g^3-g'(0)^2=0$, where we have used $f'(0)=0$. Next we multiply the equation by $(-\frac{a}{6})^2$ and substitute $w=-\frac{a}{6}g$ and we get ${w'}^2-4w^3-w'(0)^2=0$. This is the canonical differential equation for the Weierstrass elliptic function $\wp$ where $g_2=0$, and $g_3=-w'(0)^2$. The general solution is $w(x)=\wp(x + x_0;0,-w(0)^2)$. Now, recalling that $f'(x)=g(x)\;$ and $f(0)=0\;$ we have: $$f(x)=-\frac{6}{a} \int^{x}_{0} \wp\big(s+x_0;\;0,-(\frac{a}{6}f''(0))^2\big) ds=\\ =-\frac{6}{a} \bigg(\zeta\big(x_0;\;0,-(-\frac{a}{6}f''(0))^2\big)-\zeta\big(x+x_0;\;0,-(-\frac{a}{6}f''(0))^2\big) \bigg)$$ where $\zeta$ is the Weierstrass zeta function. Moreover we have $$ f'(x)=-\frac{6}{a} \wp\big(x+x_0;\;0,-(-\frac{a}{6}f''(0))^2\big)$$ and $$ f''(x)=-\frac{6}{a} \wp'\big(x+x_0;\;0,-(-\frac{a}{6}f''(0))^2\big)$$

We can use the identity for $f'(x)$ with $x=0$ to determine $x_0$ recalling that the boundary conditions require $f'(0)=0$: $$x_0=\wp^{-1}\big(0;\;0, -(\frac{a}{6}f''(0))^2\big) $$ We need another relation to determine $f''(0)$ having $f(1)=1=f''(1)$. Taking definitions of $f(x)$ and $f''(x)$ we can in principle determine $f''(0)$. We get: $$1=-\frac{6}{a}\bigg( \zeta\big(x_0;\;0,-(-\frac{a}{6}f''(0))^2\big)-\zeta\big(1+x_0;\;0,-(-\frac{a}{6}f''(0))^2\big) \bigg)$$ and $$1=-\frac{6}{a}\wp'\big(1+x_0;\;0,-(-\frac{a}{6}f''(0))^2\big) $$ Mathematica cannot solve these equations symbolically, even much simpler ones, e.g.

Reduce[2 == WeierstrassPPrime[x, {0, -4}] && 0 < x < 10, x]

Reduce::nsmet: This system cannot be solved with the methods available to Reduce.

although one can find it simply with e.g. FindRoot. When working with numeric functions like FindRoot one should be very careful to avoid possible problems with appropriate domains of functions. Let's solve numerically the identity involving $\zeta$ with respect to $f''(0)$ for a parameter $a$ (with the former choice of branch $a<0$): and the solution (with all parameters symbolic but one - $g_3=-(-\frac{a}{6}f''(0))^2$)

With[{a = -7.82}, 
  w02 =(-(a/6) fb0)^2 /. Chop@
  FindRoot[-(a/6) == WeierstrassZeta[ InverseWeierstrassP[0, {0, -(a/6 fb0)^2}],
                                      {0, -(a/6 fb0)^2}] -
           WeierstrassZeta[ 1 + InverseWeierstrassP[0, {0, -(a/6 fb0)^2}],
                            {0, -(a/6 fb0)^2}],
  { fb0, 1}];
  Plot[-(a/6) (WeierstrassZeta[ InverseWeierstrassP[0, {0, -w02}],
                               {0, -w02}] - 
               WeierstrassZeta[ x + InverseWeierstrassP[0, {0, -w02}],
                               {0, -w02}]),
        {x, -5/4, 5/4}, PlotStyle -> Thick]]

This is the way to work more or less symbolically with elliptic functions.

Answer 2

Comment

You need to use the correct syntax for the built-in function DSolve.

Apart from this, DSolve is only able to give a general solution not a particular one.

Eq = f'''[x] + a*f'[x]^2 == 0

DSolve[{Eq}, f[x], x]

DSolve[{Eq, f[0] == 0, f'[0] == 0, f[1] == 1}, f[x], x]

DSolve::bvfail: For some branches of the general solution, unable to solve the conditions.

Only one condition can be utilized,

DSolve[{Eq, f[0] == 0}, f[x], x]

So, the other choice is to go for numerical solution if you are interested in, using NDSolve.

sol = NDSolve[{Eq /. {a -> 1}, f[0] == 0, f'[0] == 0, f[1] == 1}, {f[x], f'[x]}, {x, 0, 1}]

Plot[Evaluate[{f[x], f'[x]} /. sol], {x, 0, 1}]

Answer 3

As above, you can use DSolve to find one constant but use Findroot to get the other constants, for a given a.

de = f'''[x] + a f'[x]^2 == 0;

DSolve[{de, f[0] == 0}, f[x], x];

% /. {C[1] -> c1, C[2] -> c2};

f[x_] = (f[x] /. %[[1]]);

g[x_] = f'[x];

For example with a = 2;

FindRoot[{f[1] - 1, g[0]} /. a -> 2, {{c1, -5}, {c2, -5}}] // Chop

{c1->-3.8528,c2->-5.35369}

ff[x_] = f[x] /. % /. a -> 2

-(-1)^(1/3) 3^(2/3)*
(-WeierstrassZeta[(-(1/3))^(1/3) (x-3.8528), {0,-5.35369}]-(0.9053-1.56803 I))

{ff[0], ff[1], ff'[x] /. x -> 0} // Chop

{0,1.,0}

If you really want to satisfy your fourth condition f''[1]==1, you can do so by also solving for a. You would then get for the constants:

{c1->-5.238,c2->-4.50964,a->0.867194}