Solving this challenging ODE

Consider the ODE:

$$w^{(4)}(x) + (L-x)w''(x) - w'(x) = 0$$

with some of the following boundary conditions:

• free: $w'' = 0$, $w'''=0$,
• clamped: $w = 0$, $w'=0$,
• pivot: $w = 0$, $w''=0$.

Two such pair of BC must be chosen: one for $x=0$ and one for $x=L$.

The goal is to find, provided 4 boundary conditions, the smallest length $L$ for which the ODE has some non-trivial solutions.

Question How to solve this problem nicely? By nicely, I would mean ideally something like NDSolve[{eqs, bcs}, w, x].

Issue with NDSolve One approach would be to solve the ODE with the BCs for a given L:

sol[L_] := NDSolve[Flatten[{EI*w''''[x] + q*(L - x)*w''[x] - q*w'[x] == 0,
w[0] == 0, w''[0] == 0, w[L] == 0, w''[L] == 0}], w, {x, 0, L}]


Unfortunately, NDSolve returns only the trivial solution (and that's not surprising).

My solution 1

I managed to find the solution for the different boundary conditions, but I'm not satisfied, because it is not very robust. It's not as direct as I would have hoped and the procedure depends on the BCs. I'll illustrate my solution for the fixed-free case ($w(0)=0$, $w'(0)=0$, $w''(L)=0$, $w'''(L)=0$) for which I use the integrated form of the governing equation.

{EI, q} = {1, 1};
eq = EI*w'''[x] + q*(L - x)*w'[x] == 0
wsol[x_] = w[x] /. First@DSolve[eq, w[x], x]
wsol2[x_] = wsol[x] /. First@Solve[wsol[0] == 0]
wsol3[x_] = wsol2[x] /. First@Solve[Limit[wsol2''[x], x -> L] == 0]
w[x_, L_] = wsol3[x]/C[1] // Simplify
toRoot[L_] = D[w[x, L], x] /. x -> 0
lengthCF = L /. FindRoot[toRoot[L], {L, 2}] // Chop
(* 1.98635 *)
Plot[w[x, lengthCF], {x, 0, lengthCF}]


Let's check if $w$ is a solution:

 zeroQ[x_] = EI*D[w[x, L], {x, 3}] + q*(L - x)*D[w[x, L], x] /. L -> lengthCF;
FindMaximum[{Abs@zeroQ[x], 0 < x < lengthCF}, {x, 0.5*lengthCF}]
(* {9.5451*10^-16, {x -> 0.993176}} *)


My solution 2

Since the ODE is linear, its solution is a linear combination of functions whose coefficients are determined by the BCs. So the problem can be solved as an eigenvalue problem, as mentioned by Michael Seifert. Since Mathematica can't solve the fourth order ODE, I consider a free BC at $x=L$ and a clamped BC at $x=0$:

eq = w'''[x] + (L - x)*w'[x] == 0;
sol = DSolveValue[eq, w, x];
bcs = {sol[0], sol'[0], Limit[1.*sol''[x] // Simplify, x -> L]}


The third BC requires to use numerical values because otherwise Mathematica wrongly believes it is indeterminate. Then it suffices to collect the coefficients to build the matrix of the eigenproblem, and find values for $L$ such that there are no trivial solutions:

mat = Coefficient[#, {C[1], C[2], C[3]}] & /@ bcs;
toRoot[L_] = Det[mat] // Simplify;
L /. FindRoot[toRoot[L], {L, 2}]
(* 1.98635 *)


That gives the same answer as previously. However, it works only for the 3d order ODE, and requires a trick for the third BC.

• Is eq missing a term? – Chip Hurst Mar 1 '18 at 23:13
• @ChipHurst No, if you take the derivative of eq you'll get the 4th order ODE. I used the integrated form because it appeared to work better (that's the kind of fiddling I'd like to avoid). – anderstood Mar 1 '18 at 23:38
• Ah, of course. Thanks – Chip Hurst Mar 1 '18 at 23:40
• I have a vague notion that you can treat this as an eigenvalue problem, which Mathematica can perhaps handle. I'll play around with it a bit and see what I can come up with. – Michael Seifert Mar 2 '18 at 20:36
• @MichaelSeifert That's right, I added a "solution 2" to illustrate this point. But I guess you were thinking of using something like DEigensystem? – anderstood Mar 2 '18 at 21:27

I have a package for solving eigenvalue boundary value problems using the Compound Matrix Method with the Evans function, which is ideally suited to this. The package is available on my GitHub, a partial explanation is at my previous question, and a good introduction to the method is in this pdf.

You can install the package using:

    PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]


Needs["CompoundMatrixMethod"]


And then we define the differential equation and some boundary conditions.

ode = w''''[x] + (L - x) w''[x] - w'[x] == 0;
freeBCs  = {w''[x], w'''[x]};
fixedBCs = {w[x], w'[x]};
pivotBCs = {w[x], w''[x]};


Fixed, Free:

This method requires that the original equation is recast into a matrix system, $$\mathbf{y}' = A \mathbf{y}$$, but my package includes a function to do that:

sys1 = ToMatrixSystem[ode, Join[fixedBCs/. x -> 0, freeBCs /. x -> L],
w, {x, 0, L}, L]


We now can evaluate the Evans function computed by the Compound Matrix Method. Zeroes of this Evans function correspond to eigenvalues of the original boundary value problem:

Plot[Evans[L, sys1], {L, 0, 5}]


and FindRoot gives the eigenvalue you found above.

FindRoot[Evans[L, sys1], {L, 2}]
(* {L -> 1.98635} *)


The same procedure will work for the other cases:

Pivot, Free:

sys2 = ToMatrixSystem[ode, Join[pivotBCs /. x -> 0, freeBCs /. x -> L],
w, {x, 0, L}, L];
FindRoot[Evans[L, sys2], {L, 3}]
(* {L -> 2.94869} *)


Free, Free:

I don't find a solution for the Free, Free case, which doesn't surprise me particularly.

Fixed, Fixed:

sys4 = ToMatrixSystem[ode, Join[fixedBCs/. x -> 0, fixedBCs /. x -> L],
w, {x, 0, L}, L];
FindRoot[Evans[L, sys4], {L, 4}]
(* {L -> 4.21019} *)


Using the function findAllRoots created by Jens in this answer, you can find all the roots with $$L \leq 20$$ with one line:

findAllRoots[Evans[L, sys4], {L, 0, 20}]
(* {4.21019, 5.39507, 6.87896, 7.90507, 9.10351, 10.0319, 11.0829, 11.943, 12.8987,
13.7072, 14.5944, 15.3621, 16.1967, 16.9307, 17.7231, 18.4289, 19.1864, 19.8678} *)


These agree with those found by MichaelE2 in his answer (I had initially missed the first one with the FindRoot above).

For completeness, I get 2.64806 for the pivot-pivot case (which agrees with Michael Seifert's solution) , and 3.74445 for the fixed-pivot case.

Note that as the original equation is non-symmetric, swapping round the end conditions gives different answers, with 3.10756 for the pivot-fixed case. The free-fixed and free-pivot cases have no solutions.

• The $qw'(x)$ term corresponds to a axial load density, like the gravity of a column. So it makes sense that the ODE is not symmetric (a column can buckle due to gravity, but not a suspended rod). It looks like a nice package! – anderstood Mar 8 '18 at 14:11
• Ok, that makes sense. I'm open for suggestions on the package if you have any when you try and use it, and my email address is on the package readme. On my to-do list is to write a paper explaining the method properly. – KraZug Mar 8 '18 at 16:24
• For those interested, this looks like a relevant reference: maths.gla.ac.uk/~xl/FYB-background.pdf. I did not read it nor try the package yet. – anderstood Mar 9 '18 at 14:38
• Glad to see confirmation of my result form the pivot-pivot case; I wasn't 100% sure that I had done everything correctly. – Michael Seifert Mar 9 '18 at 20:45
• For the fixed/fixed case, george2079 and I got a smaller value for L than 5.39507. (I got 5.39507 as the next smallest.) Might be a starting value problem for FindRoot. – Michael E2 Mar 13 '18 at 15:17

Chebyshev Series solution

A change of variables $\left\{ x \mapsto \frac{L}{2} (t+1)\,,\ w(x) \mapsto \frac{L}{2} u(t)\right\}$ converts the OP's differential equation to

$$(t-1) u''(t)+u'(t)=\left(\frac{2}{L}\right)^3 u^{(4)}(t),\quad -1\le t\le 1$$

Integrating four times, we get (with $\int^k$ denoting repeated integration with respect to $t$) the eigenvalue problem:

$${\cal L}u(t) \equiv-3\int^3 u(t) + (t-1) \int^2 u(t) =\left(\frac{2}{L}\right)^3 u(t)$$

Let $$u(t) = \sum_j^\infty c_j\,T_j(t),\quad -1 \le t \le 1,$$ be the Chebyshev series expansion of $u$. For a function $u$ that is analytic in a neighborhood of $-1 \le t \le 1$, the coefficients $c_j$ eventually decrease toward zero geometrically. Therefore it is possible to achieve as accurate an approximation as desired by truncating the series at a high enough degree $N$. The error is bounded by $$\sum_{j=N+1}^\infty |c_j| \approx |c_{N+1}|\big/(1-\rho)\,,$$ where $\rho$ is the geometric rate of convergence. We just need to compute the Chebyshev series of ${\cal L}u(t)$, assuming the coefficients for degree greater than $N$ are zero. This will yield an $N+1$ dimensional linear operator whose real positive eigenvalues and corresponding eigenvectors represent solutions to the OP's problem.

The basic workflow

There are some Chebyshev utilities given at the end of this answer:

• chebFunc[c, {a, b}] represents a function with Chebyshev coeffiecients c = {c0, c1,..., cN}, over the interval {a, b}; y = chebFunc[c, {a, b}][x] evaluates the function.
• chebPlus[a, b] adds two Chebyshev series.
• chebTimes[a, b] multiplies two Chebyshev series.
• chebDerivative[c, {a, b}] differentiates the Chebyshev series c scaled over the interval {a, b}.
• chebIntegrate[c, {a, b}, k] integrates the Chebyshev series c, plus k.

The ODE

Block[{EI, q},
{EI, q} = {1, 1};
eq = EI*w''''[x] + q*(L - x)*w''[x] - q*w'[x] == 0
]
(*  -Derivative[1][w][x] + (L - x)*Derivative[2][w][x] + Derivative[4][w][x] == 0  *)

eq /. {w -> (L/2 u[(2/L # - 1)] &), x -> L (t + 1)/2} // ExpandAll
(*  -u'[t] + u''[t] - t u''[t] + (8 u''''[t])/L^3 == 0  *)


I did this check in my head originally, since Integrate puts the brakes on integration by parts when one of the factors becomes an unevaluated integral.

Nest[Integrate[# // Expand, t] //.
{HoldPattern[Integrate[t f_, t]] :>
Identity[t Integrate[f, t] - Integrate[Integrate[f, t], t]],
HoldPattern[Integrate[p_Plus, t]] :> (Integrate[#, t] & /@ p)} &,
-(-u'[t] + u''[t] - t u''[t]), 4] // Simplify


This and eq above are important:

(* convert u solution to w solution *)
wsol[x_, L_, u_: usol] := L/2 u[(2/L*x - 1)]


Worked example

Set up the Chebyshev series cj, its derivatives and integrals. The derivatives are used for the boundary conditions and the integrals are used to set up the linear integral operator ${\cal L}$, the basis for which I called ode for some reason. The coefficient array of the linear component of ode is the matrix for ${\cal L}$, whose eigenvalues and eigenvectors solve the problem. The difference in the first few eigenvalues for degrees 32 and 64 differ very little. The differences can be used to determine convergence. The smaller ones vary quite a bit, and are either not very close to their true values or are side effects of truncating the series.

The boundary conditions translate to linear relations on the Chebyshev coefficients. Each can be used to eliminate a Chebyshev coefficient. Which ones varies by types of BCs. This leads to the need for getCoeffForBCs[], which was determined by inspection.

Clear[cjD, cjI];
degree = 64;
cj = Table[c[j], {j, 0, degree}];  (* Chebyshev series (coefficients) *)
dcj = Rest@NestList[chebDerivative[#, {-1, 1}] &, cj, 3];
icj = Rest@NestList[Expand@chebIntegrate[#, {-1, 1}, 0] &, cj, 4];

(* derivatives (at boundaries) and integrals (indefinite) *)
cjD[0, t : -1 | 1] := cj.t^Range[0, Length@cj - 1];
cjD[n_Integer?Positive, t : -1 | 1] := dcj[[n]].t^Range[0, Length@dcj[[n]] - 1];
cjI[n_Integer?Positive] := icj[[n]];

(* BCs in Chebyshev series *)
Clear[freeBC, clampedBC, pivotBC];
freeBC[t : -1 | 1] :=    {cjD[2, t], cjD[3, t]};  (* derivatives 2, 3 *)
clampedBC[t : -1 | 1] := {cjD[0, t], cjD[1, t]};  (* derivatives 0, 1 *)
pivotBC[t : -1 | 1] :=   {cjD[0, t], cjD[2, t]};  (* derivatives 0, 2 *)

(* code depends on global variable cj *)
ClearAll[getBCs, getCoeffForBCs];
SetAttributes[getCoeffForBCs, Orderless];
getCoeffForBCs[freeBC, pivotBC] := {1, 3, 4, 5};
getCoeffForBCs[freeBC, freeBC] := {3, 4, 5, 6};
getCoeffForBCs[_, _] := {1, 2, 3, 4};
getBCs[{left_, right_}] := First@Solve[Join[left[-1], right[1]] == 0,
cj ~Part~ getCoeffForBCs[left, right]]

(* Here's where we solve a particular problem *)
bcs0 = {clampedBC, clampedBC};  (* choose type {left BC, right BC} *)
bcs = getBCs@bcs0;

(* obtained by integrating four times *)
ode = chebPlus[chebTimes[-{1, -1}, cjI[2] /. bcs], -3 cjI[3] /. bcs];

ca = CoefficientArrays[ode, cj]  (* ca[[2]] is the desired linear operator *)


Check the matrix. The BCs determine & eliminate four of the coefficients. The coefficients beyond degree will be treated as zero, and we'll drop the last few rows and square up the matrix.

MatrixPlot@ca[[2]]


Get eigenvalues and select positive real ones. Some of the smaller ones tend to change with the degree; one might choose a lower limit such as 0.001 below.

evv = Eigenvalues[N[ca[[2]][[1 ;; Length@cj, 1 ;; Length@cj]], MachinePrecision]] //
Chop[#, 10^-40] &;  (* nix the zero imaginary parts *)

Cases[evv, v_Real /; v > 0.001]
(*
{0.107198, 0.0509448, 0.0245766, 0.0161947, 0.0106039, 0.00792391,
0.00587671, 0.00469627, 0.00372781, 0.00310627, 0.00257353, 0.00220668,
0.00188284, 0.0016484, 0.00143704, 0.00127818, 0.00113269, 0.00102009}
*)

Quiet[Solve[(L/2)^3 == 1/#, Reals], {Solve::nddc, Power::infy, Solve::naqs}] & /@
Cases[evv, v_Real /; v > 0.001];
Flatten[DeleteCases[%, _Solve], 1]
(*
{{L -> 4.21019}, {L -> 5.39507}, {L -> 6.87896}, {L -> 7.90507}, {L -> 9.10351},
{L -> 10.0319}, {L -> 11.0829}, {L -> 11.943}, {L -> 12.8987}, {L -> 13.7072},
{L -> 14.5944}, {L -> 15.3621}, {L -> 16.1967}, {L -> 16.9307}, {L -> 17.7231},
{L -> 18.4289}, {L -> 19.1864}, {L -> 19.8678}}
*)


Note I got @kraZug's 5.39507 but I also got a smaller one 4.21019, for which I could compute the solution with NDSolve. (I just noticed @george2079 got the smaller one, too.) For the rest of the problems, I get the same as @kraZug:

Constructing a solution from an eigenvector

We take the coefficient vector cj and substitute the boundary condition relations bcs and the components of the eigevector:

(* the solution's Chebyshev series is given by the eigenvector *)
usol = chebFunc[
cj /. bcs /.
Chop[Last@ Eigenvectors[N[ca[[2]][[1 ;; Length@cj, 1 ;; Length@cj]], 32], 1],
10^-44]],
{-1, 1}];


This gives the solutions in terms of $u(t)$. To get $w(x)$, we can call it as follows:

wsol[x, L, usol]


We can use it to check lower eigenvalue 4.21019 in the clamped/clamped BVP. The boundary values are sensitive in the way the roots of $(x-a)^2$ are sensitive to perturbations $(x - a)^2 \pm \epsilon$: You can only expect at best an accuracy of about half working precision. FindRoot expects 10^-8; to avoid a warning message, I lowered AccuracyGoal from 8 to 7. (One could also increase WorkingPrecision, but you have to raise it throughout the code below.)

(* store the eigenvalue in L0 *)
With[{ev = Last@Eigenvalues[N[ca[[2]][[1 ;; Length@cj, 1 ;; Length@cj]]], 1]},
L0 = L /. First@Solve[(L/2)^3 == 1/ev, Reals]
];


If L0 = 4.21019 is an eigenvalue, we can fix w''[0] at any nonzero value and shoot for the BCs at x = L0.

(* shooting for the BVP *)
Block[{obj, L},
L = L0;
obj[p2_?NumericQ, p3_?NumericQ] :=
Norm@{w[L], w'[L]} /.
(wCC = NDSolve[Flatten[{eq, w[0] == 0, w'[0] == 0, w''[0] == p2, w'''[0] == p3 }],
w, {x, 0, L}]);
w3 = FindRoot[obj[10, p3], {p3, 3}, AccuracyGoal -> 7, PrecisionGoal -> 8];
{w[L], w'[L]} /. wCC
]
(*  {{2.76458*10^-8, -6.66278*10^-8}}  *)

With[{scale = (w[1] /. First[wCC])/wsol[1, L0, usol]},
Show[
ListLinePlot[w /. First[wCC], PlotStyle -> AbsoluteThickness[3]],
Plot[scale * wsol[x, L0, usol], {x, 0, L0},
PlotLabel -> Row[{"L = ", N@L0, " ", bcs0}], PlotStyle -> {White, Dashed}]
]
]


Routine for OP's problems

(* uses global variables cj, dcj, icj *)
Clear[minLength, setupODE, allEigs];
setupODE[degree_] := (
cj = Table[\[FormalC][j], {j, 0, degree}];  (* Chebyshev series (coefficients) *)
dcj = Rest@NestList[chebDerivative[#, {-1, 1}] &, cj, 3];
icj = Rest@NestList[Expand@chebIntegrate[#, {-1, 1}, 0] &, cj, 4]
);
allEigs[{left_, right_}, n_] := Module[{bcs, ode, ca},
setupODE[n];
bcs = getBCs[{left, right}];
ode = chebPlus[chebTimes[-{1, -1}, cjI[2] /. bcs], -3 cjI[3] /. bcs];
ca = CoefficientArrays[ode, cj];
Eigenvalues[N@ca[[2]][[1 ;; n + 1, 1 ;; n + 1]]] //
Chop[#, 10^-40] &
];

(* returns smallest length L *)
minLength[bcTypes : {left_, right_}, n_] :=
Module[{bcs, ode, ca, ev, res},
ev = Max@Cases[allEigs[bcTypes, n], v_Real /; v > 0.001]; (* spurious(?) small eigs *)
If[ev > 0, (* N.B. Max[{}] = -Infinity *)
res = Solve[(L/2)^3 == 1/ev, Reals];
First@res,
"No pos. eigenvalues"]
];

(* returns all lengths L - longer ones may be spurious *)
ClearAll[allLengths];
allLengths[bcTypes : {left_, right_}, n_] :=
Module[{bcs, ode, ca, ev, res},
ev = Cases[allEigs[bcTypes, n], v_Real /; v > 0];
If[Length@ev > 0,
res = Solve[(L/2)^3 == 1/#, Reals] & /@ ev;
Flatten[res, 1],
"No positive eigenvalues"]
];


The results in the table above may be obtained with the following:

bvpLengths = AssociationMap[minLength[#, 64] &,
{{freeBC, freeBC},    {freeBC, clampedBC},    {freeBC, pivotBC},
{clampedBC, freeBC}, {clampedBC, clampedBC}, {clampedBC, pivotBC},
{pivotBC, freeBC},    {pivotBC, clampedBC},  {pivotBC, pivotBC}}];


Chebyshev utilities

These operations on Chebyshev series are based on standard properties of Chebyshev polynomials.

(*
* Chebyshev series utilities
*)
ClearAll[chebFunc];
chebFunc::usage =
"f = chebFunc[c,{a,b}], c = {c0, c1,..., cn} Chebyshev coefficients, over the interval {a,b}; y = chebFunc[c,{a,b}][x] evaluates the function";
chebFunc[c_, {a_, b_}][x_] := chebFunc[c, {a, b}, x];
chebFunc[c_?(VectorQ[#, NumericQ] &), {a_?NumericQ, b_?NumericQ}, x_?NumericQ] :=
ChebyshevT[Range[0, Length[c] - 1], (2 x - (a + b))/(b - a)].c;
chebFunc[c_?(VectorQ[#, NumericQ] &), {a_?NumericQ, b_?NumericQ}, x_?(VectorQ[#, NumericQ] &)] :=
Cos[Outer[Times, ArcCos[(2 x - (a + b))/(b - a)], Range[0, Length[c] - 1]]].c;
chebFunc /: Normal[chebFunc[c_?VectorQ, {a_, b_}, x_]] :=
Evaluate@ChebyshevT[#, (2 x - (a + b))/(b - a)] &[Range[0, Length[c] - 1]].c;

(* arithmetic and calculus on Chebyshev series *)

ClearAll[chebPlus, chebTimes, chebDerivative, chebIntegrate];
chebPlus::usage = "chebPlus[a, b] adds two Chebyshev series.";
chebTimes::usage =
"chebTimes[a, b] multiplies two Chebyshev series.";
chebPlus[aa_?VectorQ, bb_?VectorQ] :=
Module[{cc},
With[{zero = If[Precision[{aa, bb}] === MachinePrecision, 0., 0]},
cc = ConstantArray[zero, Max[Length[aa], Length[bb]]]
];
cc[[;; Length[aa]]] = aa;
cc[[;; Length[bb]]] += bb;
cc
];
chebTimes[aa_?VectorQ, bb_?VectorQ] := Module[{cc},
With[{zero = If[Precision[{aa, bb}] === MachinePrecision, 0., 0]},
cc = ConstantArray[zero, Length[aa] + Length[bb] - 1]
];
Do[
With[{dc = aa[[i]] bb},
cc[[i ;; i + Length@bb - 1]] += dc;
If[i < Length@bb,
cc[[1 ;; i]] += dc[[i ;; 1 ;; -1]];
cc[[2 ;; 1 + Length@bb - i]] += dc[[i + 1 ;; Length@bb]],
cc[[1 - Length@bb + i ;; i]] += Reverse@dc;
]
],
{i, Length@aa}];
cc/2
];
chebDerivative::usage =
"chebDerivative[c, {a,b}] differentiates the Chebyshev series c scaled over the interval {a,b}";
chebDerivative[c_, {a_, b_}] := Module[{c1 = 0, c2 = 0, c3},
2/(b - a) MapAt[#/2 &,
Reverse@ Table[
c3 = c2; c2 = c1;
c1 = 2 (n + 1)*c[[n + 2]] + c3, {n, Length[c] - 2, 0, -1}],
1]
];
chebIntegrate::usage =
"chebIntegrate[c, {a, b}, k] integrates the Chebyshev series c, plus k";
chebIntegrate[c0_, {a_, b_}, k_: 0] := Module[{c, i, i0},
c[1] = 2*First[c0];
c[n_] /; 1 < n <= Length[c0] := c0[[n]];
c[_] := 0;
i = (b - a)/2 Table[(c[n - 1] - c[n + 1])/(2*(n - 1)), {n, 2, Length[c0] + 1}];
i0 = i[[2 ;; ;; 2]];
Prepend[i, k - Sum[(-1)^n*i0[[n]], {n, Length[i0]}]]
];

• Chebyshev methods are awesome. Maybe you should start collecting these into a package. – J. M.'s technical difficulties Mar 13 '18 at 15:26
• @J.M. Well, I have started.... – Michael E2 Mar 18 '18 at 15:17
• On that note: apart from Boyd, have you already seen Trefethen's book? There's so much in there I'd like to try doing in Mathematica, but not enough time... – J. M.'s technical difficulties Mar 18 '18 at 15:25
• @J.M. Yes. I haven't done much of the later stuff. I implemented Carathéodory-Fejér, which did fantastic on Sin[x] over {0, 10}, pretty much converting the Chebyshev approximation to the minimax one. But it made the approximation to the implicit function x == y Exp[10 y^2] over {0, 2} worse; I've yet to find out why. – Michael E2 Mar 18 '18 at 15:49

Here is how to get the result directly from @Vsevolod analytic solution:

seeking a solution to the case w[0]==w'[0]==w''[L]==w'''[L]==0, substituting w'[x]=v[x] :

vf[x_, L_, c1_, c2_] =
DSolveValue[{v'''[x] + (L - x) v'[x] - v[x] == 0, v[0] == 0}, v[x],
x] /. {C[1] -> c1, C[2] -> c2};


now for w'''[L] == v''[L] ==0:

 Limit[(D[vf[x, L, c1, c2], {x, 2}]), {x -> L}]


we get an expression involving only c1 which then must be zero. (notice v''[L] == 0 could have been given to DSolve as well).

now for w''[L] == v'[L] ==0 :

  D[vf[x, L, 0, c2], x] /. x -> L // Simplify


-((3^(1/6) c2 AiryAi[-L] + (c2 AiryBi[-L])/3^(1/3))/( AiryBi[-L] Gamma[1/3]))

c2==0 gives the trivial solution, so we look for L to make this zero.

el = L /. FindRoot[%/c2 == 0, {L, 2}][[1]]


1.98635

this simplifies (by hand) to FindRoot[BesselJ[-1/3,2 L^(3/2)/3], {L, 2}]

the resulting v is linear in c2

vf[x, el, 0, c2] // Chop // Simplify


c2 (AiryAi[-1.9863527074304728+ x] + 0.5773502691896263 AiryBi[-1.9863527074304728 + x])

which is (by hand)

 2/3 c2 Sqrt[L - x] BesselJ[-1/3, 2/3 (L - x)^(3/2)]


c2 can be normalized out and integrate to get w:

NDSolveValue[{w'[x] == vf[x, el, 0, c2]/c2, w[0] == 0}, w[x], {x, 0, el}]
Plot[%, {x, 0, el}]


note we can use FindRoot to find higher order solutions:

L=5.29562

edit: the w[0]==w'[0]==w[L]==w'[L]==0 problem:

vf[x_, L_, c1_, c2_] =
DSolveValue[{v'''[x] + (L - x) v'[x] - v[x] == 0, v[0] == 0,
v[L] == 0}, v[x], x] /. C[2] -> c2 // Simplify;


to have w[0]==w[L] the integral must be zero:

g[L_?NumericQ] := NIntegrate[vf[x, L, c1, c2]/c2, {x, 0, L}]
el = L /. FindRoot[g[L] == 0, {L, 4.2}]


4.21019

NDSolveValue[{w'[x] == vf[x, el, 0, c2]/c2, w[0] == 0},
w[x], {x, 0, el}]
Plot[%, {x, 0, el}, PlotRange -> All]


first 3 solutions..

• Instead of FindRoot, you can also use NSolve[BesselJ[-1/3, 2 L^(3/2)/3] == 0 && 0 < L < 10, {L}]. (+1) – Michael E2 Mar 13 '18 at 15:08
• As an alternative to @Michael's suggestion: N[(3 BesselJZero[-1/3, Range[7]]/2)^(2/3)] – J. M.'s technical difficulties Mar 14 '18 at 1:46

NDEigenvalue method (partial answer)

With a bit of massaging, the pivot-pivot case can be treated in Mathematica using NDEigenvalue. We have to start by rewriting the equation in a quasi-eigenvalue form. Defining $v = w''$, we can replace the fourth-order equation by two second-order equations: \begin{align*} -v'' + x v + w' &= L v \\ w'' - v &= 0 \end{align*} Now, NDEigenvalue expects to solve for a system of equations of the form $$\mathcal{L}_1 (v, w) = L v \\ \mathcal{L}_2 (v,w) = L w$$ which is not quite what we have here. But we can kludge this by trying to solve the system \begin{align*} -v'' + x v + w' &= L v \\ M(w'' - v) &= L w \end{align*} If $M \gg L$, the solutions to these equations will approximately satisfy the original equations we have. Our strategy will therefore be to pick a large number $M$ and use NDEigenvalue to solve the eigenvalue problem; if the result is that $L \ll M$, then we have a good approximate to the solution of our original problem.

The other difficulty is dealing with the boundary conditions, which also depend on $L$. We address this by solving the eigenvalue problem on the region $[0, a]$ instead, for general $a$. The critical value of $a$ will be the value for which the lowest eigenvalue $L$ of the system is equal to $a$. We can use numerical root-finding techniques to find this critical value.

The easiest case to deal with, given the default usage of NDEigenvalue, is the pivot-pivot case. In this case, the boundary conditions for our second-order system are simple Dirichlet boundary conditions: $$w(0) = v(0) = w(a) = v(a) = 0.$$

Code:

m = 10^8;
eqns = {-v''[x] + x v[x] + w'[x], m (w''[x] - v[x])};
lowesteigenval[a_] :=
NDEigenvalues[
Join[eqns, {DirichletCondition[v[x] == 0, True],
DirichletCondition[w[x] == 0, True]}], {v[x], w[x]}, {x, 0, a}, 1]


Playing around with the lowesteigenval function shows us that the critical value of a is between about 2 and 3, which allows us to call FindRoot with the secant method:

critval = FindRoot[lowesteigenval[L] - L, {L, 2, 3}]
(* {L -> 2.64806} *)


We can then use this critical value to plot the solution for $w(x)$. (Note that the way we have defined things, $w(x)$ is the second function returned by NDEigensystem.

soln = NDEigensystem[
Join[eqns, {DirichletCondition[v[x] == 0, True],
DirichletCondition[w[x] == 0, True]}], {v[x], w[x]}, {x, 0, L /. critval}, 1]
Plot[soln[[2, 1, 2]], {x, 0, L /. critval}]


As a sanity check, increasing or decreasing the value of m does not substantially change this answer (until one starts bumping up against machine precision issues.)

Further work

The free and fixed boundary conditions are probably accessible by this method as well. However, under this method, the boundary conditions are not pure Dirichlet BCs. for a clamped end we have $w(a) = w'(a) = 0$, while for a free end we have $v(a) = v'(a) = 0$. In other words, we have a mix of Dirichlet and Neumann boundary conditions, and NDEigensystem and its related functions handle Neumann BCs rather differently from Dirichlet BCs in their syntax. I suspect that it's possible to extend this method to the other types of boundary conditions; if I have further time, I'll return to this answer and see if I can extend it.

Your equation can be solved analytically.

$w^{\{4\}}+(L−x)w′′−w′=0$

$w'=v$

$v^{\{3\}}+(L−x)v′−v=0$

Note that

$(L-x)v'-v=((L-x)\cdot v)'$

Thus

$v^{\{3\}}+((L-x)\cdot v)'=0$

$v^{\{2\}}+((L-x)\cdot v)=C$

The later one is processed by mathematica.

DSolve[v''[x] + (L - x) v[x] == M, v[x], x]

• It can be solved analytically for $w'$, but not for $w$. Also, Solve[v'[L] == 0 && v'[0] == 0] can't be solved for example. How would you find the critical $L$? – anderstood Mar 4 '18 at 12:30
• @anderstood Yes it can DSolve[{v''[x] + (L - x) v[x] == M, v'[L] == 0, v'[0] == 0}, v[x], x] – Vsevolod A. Mar 4 '18 at 13:08
• DSolve also can solve the third-order ODE. – bbgodfrey May 17 '18 at 0:09

Fourier methods

I will only discuss the pivot-pivot case since bases for this cases are easiest to determine. The other cases will "just" need "a bit" more linear algebra. To this end, I will cast the problem into a linear, generalized eigenvalue problem.

Step 1 (is the same for all boundary conditions)

The left hand side of the equation looks like this:

LHS = (D[w[x], {x, 4}] + (L - x) D[w[x], {x, 2}] - D[w[x], x]);


Via substitutions

subs = x == a t + b /.
Solve[{0 == a t + b /. {t -> -Pi}, L == a t + b /. {t -> Pi}}, {a, b}][[1]];
t2x = Solve[subs, t][[1]]
x2t = Solve[subs, x][[1]]
w[x_] := Evaluate[u[t] /. t2x]
Lsubs = Solve[L/(2 Pi) == Power[λ, 1/3], L][[1]];
LHS2 = ((L/(2 Pi))^4 LHS /. x2t // Simplify // Together) /. Lsubs //
Expand


$\left\{t\to -\frac{\pi (L-2 x)}{L}\right\}$

$\left\{x\to \frac{L (t+\pi )}{2 \pi }\right\}$

$\left\{L\to 2 \pi \sqrt[3]{\lambda }\right\}$

we arrive at this:

w[x_] := Evaluate[
u[t] /. t2x]; LHS2 = ((L/(2 Pi))^4 LHS /. x2t // Simplify // Together) /. Lsubs // Expand


$u^{(4)}(t) + \lambda \, \big( (\pi -t) \, u''(t) - u'(t) \big)$

With the differential operators

$$A_4 = \Big(\frac{\operatorname{d}}{\operatorname{d} t} \Big)^4, \quad A_2 = (t - \pi) \, \Big(\frac{\operatorname{d}}{\operatorname{d} t} \Big)^2, \quad \text{and} \quad A_1 = \frac{\operatorname{d}}{\operatorname{d} t},$$

we are essentially seeking for eigenvalues $\lambda$ fulfilling

$$A_4 \, u - \lambda \, (A_2 + A_1) \, u = 0$$

for a function $u \colon [-\pi , \pi ] \to \mathbb{R}$ with

$$u(-\pi) = u(\pi) = 0 \quad \text{and} \quad u''(-\pi) = u''(\pi) = 0.$$

Discretization

Fine, the functions $e_k(t) = \sin(k(t+\pi)/2)$ form a complete $L^2$-orthonormal system. The matrix entries

$$(A_1)_{k,m} = \int_{-\pi}^\pi \sin(k(t+\pi)/2) \dfrac{\mathrm{d}}{\mathrm{d}t} \sin(m(t+\pi)/2) \, \operatorname{d} t$$

$$(A_2)_{k,m} = \int_{-\pi}^\pi \sin(k(t+\pi)/2)(t-\pi) \dfrac{\mathrm{d}^2}{\mathrm{d}t^2} \sin(m(t+\pi)/2) \, \operatorname{d} t$$

$$(A_4)_{k,m} = \int_{-\pi}^\pi \sin(k(t+\pi)/2) \dfrac{\mathrm{d}^4}{\mathrm{d}t^4} \sin(m(t+\pi)/2) \, \operatorname{d} t$$

can be easily computed symbolically:

buildMatrix[coef_, order_, n_] := Block[{offdiag, diag, f},
offdiag = Simplify[Integrate[
Sin[k*(t + Pi)/2]*coef*D[Sin[m*(t + Pi)/2], {t, order}], {t, -Pi,
Pi}], k \[Element] Integers && m \[Element] Integers];
diag = Simplify[Integrate[
Sin[k*(t + Pi)/2]*coef*D[Sin[k*(t + Pi)/2], {t, order}], {t, -Pi,
Pi}], k \[Element] Integers];
f = {m, k} \[Function] Evaluate@If[m == k, Evaluate@diag, Evaluate@offdiag];
Array[f, {n, n}]]


For 10 basis functions, we obtain:

n = 10;
A1 = buildMatrix[1, 1, n];
A2 = buildMatrix[t - Pi, 2, n];
A4 = buildMatrix[1, 4, n];


Then eigenvalues are given by:

eigs = Sort[Eigenvalues[N@{Normal@A4, Normal[A2 + A1]}]];
Lsubs /. \[Lambda] -> eigs[[1]]
(* {L -> 2.64806} *)


which is the same result as in the other answers.

Remarks

For clamped boundary conditions on both side, functions like $(-1)^{k+2} k \sin (t)+\sin (k t)$ or $(-1)^{k+2} \cos (t)+\cos (k t)$ can be used and so an. Admittedly, things start to become quite confusing when different boundary conditions at the two ends are used...

• Shouldn't $A_1$ and $A_4$ be defined by $\int_{-\pi}^\pi \sin(kt)\sin^{(\alpha)}(kt)$ with $\alpha \in \{1,4\}$? Or would there be a missing factor $\pi$ in the definition of $A_1$ and $A_4$? – anderstood Mar 13 '18 at 18:47
• The diagonal entries should look like $\int_{-\pi}^\pi \sin(kt) \left( \tfrac{d}{d t}\right)^\alpha \sin(kt) \, d t$ such that differentiating is essentially equivalent to multiplication with frequencies. But thinking of it: this would turn $A_1$ into the zero matrix. Thinking about it even more: At least the $\cos$ function oughts to be also a member of the basis. Well, I have to admit: The whole attempt is pretty broken... – Henrik Schumacher Mar 13 '18 at 19:30
• The strange thing is that i) the result would not depend on $n$, ii) this would lead to a closed-form expression for $L_\mathrm{min}$, which mustn't exist (I'm not 100% sure in the pivot-pivot case). – anderstood Mar 13 '18 at 19:32
• OK I spotted the mistake: $\sin(kt)$ is not a basis, instead use $\sin(k(t+\pi)/2)$ and then it works! I got $2.64806$ with $n=500$, $2.64827$ with $n=3$. Great!! As per the $\cos$, if you include them in the basis, they will "disappear" because of the BCs, so using only sines is fine. – anderstood Mar 13 '18 at 20:39
• Huh? It works? Wow, that's great! I did't believe that it can be fixed. Hm. For k=1, I get $\sin(k(t+\pi)/2) = \cos(t/2)$. That must be the $\cos$-function that I meant (but was not able to nail it down)... Now I realize that this is essentially the basis of the discrete $\sin$-transform...Very good job! – Henrik Schumacher Mar 13 '18 at 21:20