How to solve fourth order differential equation?

Question

I have a differential equation of this type:

y[x] - 1 - 2*l^2*y''[x] + l^4*y''''[x] == 0

(where l is a parameter and l>=0),

with boundary condtions that needs to be satisfied:

y[0] == 0, y'[0] == 0, y'[±∞] == 0, y''[±∞] == 0

The solution to this differential equation is:

y[x_] := 1 - Exp[-Abs[x]/l]*(1 + Abs[x]/l)

but I can't get it right. (I added a graph of solution function for l=0.5)

When I try to use DSolve (with 4 initial conditions because it's a fourth ordered differential equation) I get an error saying "For some branches of the general solution, unable to compute the limit at the given points. Some of the solutions may be lost."

DSolve[{y[x] - 1 - 2*l^2*y''[x] + l^4*y''''[x] == 0, y[0] == 0, 
        y'[0] == 0, y'[∞] == 0, y''[∞] == 0}, y[x], x]

I don't understand how this solution was obtained.

Can anyone help me?

Answer 1

When $l=0$, the solution is clearly $y(x)=1$:

eq = y[x] - 1 - 2 l^2 y''[x] + l^4 y''''[x] == 0;  
Solve[eq /. l -> 0, y[x]]
(* {{y[x] -> 1}} *)

And it's the limit for the given solution:

Limit[1 - Exp[-Abs[x]/l]*(1 + Abs[x]/l), l -> 0, Assumptions -> x != 0]
(* 1 *)

Then let's focus on the $l>0$ case. We know DSolve cannot handle boundary condition (b.c.) at infinity directly, so we first solve with the 2 b.c.s at $x=0$:

generalsol = y[x] /. First@DSolve[{eq, y[0] == 0, y'[0] == 0}, y, x]

(* -((E^(-(x/l)) (-E^((x/l)) l + E^((2 x)/l) l - E^((2 x)/l) x - l C[1] + 
      E^((2 x)/l) l C[1] - 2 E^((2 x)/l) x C[1] - l x C[2] + E^((2 x)/l) l x C[2]))/l) *)

Calculate $y'(x)$ and transform it a bit:

Collect[D[generalsol, x], Exp[_]]
(* E^(x/l) (x/l^2 + C[1]/l + (2 x C[1])/l^2 - C[2] - (x C[2])/l) + 
   E^(-(x/l)) (-(C[1]/l) + C[2] - (x C[2])/l) *)

Let's consider $x\in(0,\infty)$ first, it's clear that the coefficient of E^(x/l) term should be 0 when $l>0$, or $y'(\infty)=0$ won't be satisfied:

cond1 = x/l^2 + C[1]/l + (2 x C[1])/l^2 - C[2] - (x C[2])/l == 0;

Similarly, we calculate $y''(x)$ and transform it a bit:

Collect[D[generalsol, x, x], Exp[_]]
(* E^(x/l) (1/l^2 + x/l^3 + (3 C[1])/l^2 + (2 x C[1])/l^3 - (2 C[2])/l - (x C[2])/l^2) + 
   E^(-(x/l)) (C[1]/l^2 - (2 C[2])/l + (x C[2])/l^2) *)

It's clear that the coefficient of E^(x/l) term should be 0, or $y''(\infty)=0$ won't be satisfied:

cond2 = 1/l^2 + x/l^3 + (3 C[1])/l^2 + (2 x C[1])/l^3 - (2 C[2])/l - (x C[2])/l^2 == 0;

So C[1] and C[2] are determined:

coef = Flatten@Solve[{cond1, cond2}, C /@ {1, 2}]
(* {C[1] -> -1, C[2] -> -(1/l)} *)

And we find the solution for $x\in(0,\infty)$:

rightsolfunc[x_] = generalsol /. coef // FullSimplify
(* 1 - (E^(-(x/l)) (l + x))/l *)

We can repeat the process above to find solution for $x\in(-\infty,0)$ of course, but here I'll skip it because it's easy to notice rightsolfunc[-x] is the solution for $x\in(-\infty,0)$, thus the complete solution is:

sol = rightsolfunc@Abs@x
(* 1 - (E^(-(Abs[x]/l)) (l + Abs[x]))/l *)