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 /.
Thread[cj ->
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]}]]
];
eq
missing a term? $\endgroup$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). $\endgroup$DEigensystem
? $\endgroup$