Discretization of 4th order ODE

I am trying to use finite difference method to find solution of the following ODE with boundary conditions : $$y^{(4)}+ P{y}^{''} = 0\;(**)$$ $$y(0) = y'(0) = 0, \;y^{''}(L) = 0,\;y^{(3)}(L)+Py'(L)=0$$ where $P$ is a positive parameter. Problem should be turned into algebraic problem for which $P_n$ should be it's eigenvalues and then use some standard methods for approximating eigenvalues. I have tried using Taylor's expansion to get approximation for the 2nd and 4th derivative : $$y{''}_i \approx\dfrac{y_{i+1}-2y_i+y_{i-1}}{h^2}$$ $$y^{(4)}_i \approx\dfrac{y_{i+2}-4y_{i+1}+6y_i-4y_{i-1}+y_{i-2}}{h^4}$$ Plugging into $(**)$: $$y_{i+2} + (-4 +Ph^2)y_{i + 1}+(6-2Ph^2)y_i+(-4+Ph^2)y_{i-1}+y_{i-2} = 0$$ I am not sure what to do with this. I should use boundary conditions and get my problem in the form $$Ay =Py$$ where A is a constant matrix.

• It's better to solve this as a system of second order equations (or even first order). Feb 16 '18 at 12:38
• I solved it already as a second order problem with the equation $-y^{''}=P(y-y(L))$. In the answer I posted below is the matrix-vector equation for second order case. Second part of the problem is to also solve it as an eigenvalue problem but with the 4th order ODE and compare results. Just for reference physical representation of the problem is a variation of en.wikipedia.org/wiki/Euler%27s_critical_load. In my case one end is fixed while the other is free to oscillate. That part is a little bit tricky since I only have boundary values at one of the endpoints.
– Tino
Feb 16 '18 at 13:07

First we fix the expression for $y_i^{(4)}$: $$y_i^{(4)}\approx\frac{y_{i+2}-4y_{i+1}+6y_i-4y_{i-1}+y_{i-2}}{h^4}$$ This doesn't work at $i=1$ because there is no $y_{-1}$ so we write out equations as best we can using the boundary conditions as much as possible: $$\begin{bmatrix}1&-1&\frac12&-\frac16&\frac1{24}&-\frac1{120}\\ 0&1&-1&\frac12&-\frac16&\frac1{24}\\ 1&0&0&0&0&0\\ 1&1&\frac12&\frac16&\frac1{24}&\frac1{120}\\ 1&2&2&\frac43&\frac23&\frac4{15}\\ 1&3&\frac92&\frac92&\frac{27}8&\frac{81}{40}\end{bmatrix}\begin{bmatrix}y_1\\hy_1^{\prime}\\h^2y_1^{\prime\prime}\\h^2y_1^{\prime\prime\prime}\\h^4y_1^{(4)}\\h^5y_1^{(5)}\end{bmatrix}=\begin{bmatrix}y_0\\hy_0^{\prime}\\y_1\\y_2\\y_3\\y_4\end{bmatrix}$$ Thus $$h^4y_1^{(4)}=-\frac{113}{12}y_0-5y_0^{\prime}+16y_1-9y_2+\frac83y_3-\frac14y_4$$ We have a similar problem for $y_{n-1}$ in that there is no $y_{n+1}$. $$\begin{bmatrix}1&1&\frac12&\frac16&\frac1{24}&\frac1{120}\\ 0&h^2P&h^2P&1+\frac12h^2P&1+\frac16h^2P&\frac12+\frac1{24}h^2P\\ 0&0&1&1&\frac12&\frac16\\ 1&0&0&0&0&0\\ 1&-1&\frac12&-\frac16&\frac1{24}&-\frac1{120}\\ 1&-2&2&-\frac43&\frac23&-\frac4{15}\end{bmatrix}\begin{bmatrix}y_{n-1}\\hy_{n-1}^{\prime}\\h^2y_{n-1}^{\prime\prime}\\h^3y_{n-1}^{\prime\prime\prime}\\h^4y_{n-1}^{(4)}\\h^5y_{n-1}^{(5)}\end{bmatrix}=\begin{bmatrix}y_n\\h^3y_n^{\prime\prime\prime}+h^2Phy^{\prime}\\h^2y_n^{\prime\prime}\\y_{n-1}\\y_{n-2}\\y_{n-3}\end{bmatrix}$$ With solution $$\left(\frac{17}4-\frac3{10}h^2P\right)h^4y_{n-1}^{(4)}=\left(4-\frac4{15}h^2P\right)y_{n-3}+\left(-15+\frac9{10}h^2P\right)y_{n-2}+18y_{n-1}+\left(3-\frac35h^2P\right)h^2y_n^{\prime\prime}+\left(-7-\frac{19}{30}h^2P\right)y_n+\left(h^3y_n^{\prime\prime\prime}+h^2Phy_n^{\prime}\right)$$ Then for $i=n$ we have $$\begin{bmatrix}1&0&0&0\\ 1&-1+\frac16h^2P&\frac1{24}&-\frac1{120}\\ 1&-2+\frac43h^2P&\frac23&-\frac4{15}\\ 1&-3+\frac92h^2P&\frac{27}8&-\frac{81}{40}\end{bmatrix}\begin{bmatrix}y_n\\hy_n^{\prime}\\h^4y_n^{(4)}\\h^5y_n^{(5)}\end{bmatrix}=\begin{bmatrix}y_n\\y_{n-1}\\y_{n-2}\\y_{n-3}\end{bmatrix}$$ And I get something like $$\left(\frac{255}4-\frac92h^2P\right)h^4y_n^{(4)}=\left(270-222h^2P\right)y_n+\left(-585+270h^2P\right)y_{n-1}+\left(360-54h^2P\right)y_{n-2}+\left(-45+6h^2P\right)y_{n-3}$$ This was all complicated enough that I almost inevitably made many mistakes so you should work it through for yourself and let me know about any discrepancies.

And I checked the first result with Excel, so maybe the expressions aren't so bad after all :)

EDIT: Now that we have done the hard stuff, we get to do the easy part. The actual equation at node $i$ is $h^4y_i^{(4)}+h^2Ph^2y_i^{\prime\prime}=0$. At node $1$ equation $1$ reads $$\left(16-2h^2P\right)y_1+\left(-9+h^2P\right)y_2+\frac83y_3-\frac14y_4=0$$ At nodes $2\le i\le n-2$ equation $i$ reads $$y_{i-2}+\left(-4+h^2P\right)y_{i-1}+\left(6-2h^2P\right)y_i+\left(-4+h^2P\right)y_{i+1}+y_{i+2}=0$$ We will promote the equation for node $n$ to be equation $n-1$: $$\left(-15+2h^2P\right)y_{n-3}+\left(120-18h^2P\right)y_{n-2}+\left(-195+90h^2P\right)y_{n-1}+\left(90-74h^2P\right)y_n=0$$ The remaining equation is quite bloody: $$\left(4-\frac4{15}h^2P\right)y_{n-3}+\left(-15+\frac{103}{20}h^2P-\frac3{10}h^4P^2\right)y_{n-2}+\left(18-\frac{17}2h^2P+\frac35h^4P^2\right)y_{n-1}+\left(-7+\frac{217}{60}h^2P-\frac3{10}h^4P^2\right)y_n=0$$ This one goes against the grain of our desired form. We will define a fictitious variable $y_{n+1}$ via equation $n$: $$\frac3{10}h^2Py_{n-2}-\frac35h^2Py_{n-1}+\frac3{10}h^2Py_n+y_{n+1}=0$$ Now we can add $h^2P$ times equation $n$ to the offending equation to get equation $n+1$: $$\left(4-\frac4{15}h^2P\right)y_{n-3}+\left(-15+\frac{103}{20}h^2P\right)y_{n-2}+\left(18-\frac{17}2h^2P\right)y_{n-1}+\left(-7+\frac{217}{60}h^2P\right)y_n+h^2Py_{n+1}=0$$ Now our system is in the form $$Ay=Bh^2Py$$ So it can be solved like an ordinary eigenvalue/eigenvector problem.

EDIT: Tacked down a sign error that was keeping everything from working. Now in Matlab it looks like

% bvp1.m

clear all;
close all;
narray = [5:100];
for k=1:length(narray),
n = narray(k);
A = zeros(n+1);
B = zeros(n+1);
A(1,1:4) = [16 -9 8/3 -1/4];
B(1,1:2) = [-2 1];
A(2,1:4) = [-4 6 -4 1];
B(2,1:3) = [1 -2 1];
for i = 3:n-2,
A(i,i-2:i+2) = [1 -4 6 -4 1];
B(i,i-1:i+1) = [1 -2 1];
end
A(n-1,n-3:n) = [-15 120 -195 90];
B(n-1,n-3:n) = [2 -18 90 -74];
A(n,n+1) = 1;
B(n,n-2:n) = [3/10 -3/5 3/10];
A(n+1,n-3:n) = [4 -15 18 -7];
B(n+1,n-3:n+1) = [-4/15 103/20 -17/2 217/60 1];
[V,D] = eig(-B\A);
[D, ind] = sort(diag(D));
P(1:3,k) = n^2*D(1:3);
end
figure;
plot(narray,P(1,:),narray,P(2,:),narray,P(3,:));
title('Lowest Eigenvalues for Buckling Beam');
xlabel('Number of Grid Points')
ylabel('\lambda^2');
legend('\lambda_1^2','\lambda_2^2','\lambda_3^2','Location','Best')
figure;
plot([1:n],V(1:n,ind(1)),'k-',[1:n],V(1:n,ind(2)),'b-',[1:n],V(1:n,ind(3)),'r-');
title('First Three Buckling Modes');
xlabel('Grid Point');
ylabel('Displacement');
legend(['\lambda_1^2 = ' num2str(P(1,end))],['\lambda_2^2 = ' num2str(P(2,end))],['\lambda_3^2 = ' num2str(P(3,end))],'Location','Best')


• This has helped tremendously!!!
– Tino
Feb 17 '18 at 12:23

I assume you were using $y_n = y(L)$ as the right boundary. Boundary conditions then simplify the expressions $$h^4y_1^{(4)}=16y_1-9y_2+\frac83y_3-\frac14y_4$$ $$\left(\frac{17}4-\frac3{10}h^2P\right)h^4y_{n-1}^{(4)}=\left(4-\frac4{15}h^2P\right)y_{n-3}+\left(-15+\frac9{10}h^2P\right)y_{n-2}+18y_{n-1}+\left(-7-\frac{19}{30}h^2P\right)y_n$$ $$\left(\frac{255}4+\frac92h^2P\right)h^4y_n^{(4)}=\left(270+222h^2P\right)y_n+\left(-585-270h^2P\right)y_{n-1}+\left(360+54h^2P\right)y_{n-2}+\left(-45-6h^2P\right)y_{n-3}$$ Last equation stays the same. For $i = 2, 3,...,n - 2$ $$y_i^{(4)}\approx\frac{y_{i+2}-4y_{i+1}+6y_i-4y_{i-1}+y_{i-2}}{h^4}$$ Combining approximations with ODE for $i =1$
$$(16-2Ph^2)y_1 + (-9+Ph^2)y_{2}+\frac{8}{3}y_3-\frac{1}{4}y_{3} = 0$$ for $i =2, 3,..., n-2$
$$y_{i+2} + (-4 +Ph^2)y_{i + 1}+(6-2Ph^2)y_i+(-4+Ph^2)y_{i-1}+y_{i-2} = 0$$ for $i = n - 1$ and $i = n$ they turn out pretty messy. Usually when I was doing problems like these of 2nd order I would either have P given and then calculate approximate solutions $y_n$, or I could make the discretization such that P would be an eigenvalue of $Ay=py$ and then get $y_n$ as the corresponding eigenvectors. For example $$2y_0-2y_1=Ph^2y_0$$ $$-y_{n+1}+2y_i-y_{i-1}=Ph^2y_i$$ $$A = \begin{bmatrix} 2 & -2 & 0 & 0 & \dots & 0 \\ 1 & 2 & -1 & 0 & \dots & 0 \\ 0 & -1 & 2 & -1 & \dots & 0\\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & 0 & -1 & 2 \end{bmatrix} \begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_n\\ \end{bmatrix}= Ph^2 \begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_n\\ \end{bmatrix}$$