# How to solve fourth order differential equation?

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, y' == 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,
y' == 0, y'[∞] == 0, y''[∞] == 0}, y[x], x]


I don't understand how this solution was obtained.

Can anyone help me? • Simpson, you've been a member for 1.5 year, please make some effort to learn formatting your post properly. Also, please show us how you tried DSolve so we can give pointed advice. – xzczd Jan 5 '19 at 12:28
• An elliptic 4th order boundary values problem with six boundary conditions sounds fishy. Btw., you can eliminate l by the substitution $x \to l \, x$. – Henrik Schumacher Jan 5 '19 at 12:48
• Also, there should be certain constraint for $l$, for example, when $l=-1$ the b.c. at infinity won't be satisfied: Limit[y'[x] /. y -> (Function[x, #1] &)[Simplify[1 - (E^(-(Abs[x]/l)) (l + Abs[x]))/l, x > 0]] /. l -> -1, x -> \[Infinity]] – xzczd Jan 5 '19 at 15:30
• @J. M . is computer-less: Could you edit y[x_] too (Abs is from the original form of the question, where 6 boundary conditions were stated.)? TIA. – user64494 Jan 5 '19 at 21:14
• @user64494 What edit are you expecting? J.M. only improved the formatting. Except for the clarification $l\geq0$, the content of this question doesn't change so far. – xzczd Jan 6 '19 at 10:35

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, y' == 0}, y, x]

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


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

Collect[D[generalsol, x], Exp[_]]
(* E^(x/l) (x/l^2 + C/l + (2 x C)/l^2 - C - (x C)/l) +
E^(-(x/l)) (-(C/l) + C - (x C)/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/l + (2 x C)/l^2 - C - (x C)/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)/l^2 + (2 x C)/l^3 - (2 C)/l - (x C)/l^2) +
E^(-(x/l)) (C/l^2 - (2 C)/l + (x C)/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)/l^2 + (2 x C)/l^3 - (2 C)/l - (x C)/l^2 == 0;


So C and C are determined:

coef = Flatten@Solve[{cond1, cond2}, C /@ {1, 2}]
(* {C -> -1, C -> -(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 *)

• Unfortunately, y[x]=1 - (E^(-(RealAbs[x]/l)) (l +Real Abs[x]))/l does not satisfy eq at $x=0$. In order to clarify my point, plot y'''[x] on $[−1,1]$. – user64494 Jan 6 '19 at 10:27
• @user64494 …Why it should? Notice this solution is just not a classical one, it's determined by 6 b.c.s while for a classical solution is just determined by 4. What's most important: this is the solution mentioned in the question i.e. the solution whose 2nd order derivative is not smooth at $x=0$ is exactly what OP is looking for. – xzczd Jan 6 '19 at 10:30
• Can you explain what you mean by "this solution is just not a classical one" and give a reference? The third derivative y'''[x] is not continuous at the origin. Therefore, the fourth derivative does not exist at the origin at all. – user64494 Jan 6 '19 at 11:17
• @user64494 For example this, see also the 1st paragraph of this page, you may search for classical solution, weak solution to find more. As to the 4th order derivative, it actually exists, though it's not smooth. Finally, once again, the non-smooth derivative at $x=0$ just doesn't matter, because this is what OP is looking for. – xzczd Jan 6 '19 at 11:58
• @user64494 …Let's forget about those terminologies for a while and go a step back. With b.c.s y == 0, y' == 0, y'[∞] == 0, y''[∞] == 0 we find the solution for $[0,\infty)$ is 1 - (E^(-(x/l)) (l + x))/l, and with y == 0, y' == 0, y'[-∞] == 0, y''[-∞] == 0 we find the solution for $(-\infty, 0)$ is 1 - (E^(x/l) (l - x))/l, then why can't we shorten the solution to 1 - (E^(-(Abs[x]/l)) (l + Abs[x]))/l? You may argue the notation is not religious, but such stuff exists everywhere. For example, how do you think about the function WhenEvent and the solutions produced by it? – xzczd Jan 6 '19 at 12:37

We may solve fourth order differential equation using NDsolve

Test input;

  l = 1;
fffp = NDSolve[{y[x] - 1 - 2*l^2*y''[x] + l^4*y''''[x] == 0,
y == 0, y' == 0, y' == 0, y'' == 0}, y[x], x];
Plot[Evaluate[y[x] /. fffp], {x, 0, 1}, PlotRange -> All] 