Discretization of 4th order ODE

Question

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.

Answer 1

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.

Wolfram alpha agrees with my last result: https://www.wolframalpha.com/input/?i=%7B%7B1,+0,+0,+0%7D,%7B1,+-1%2B1%2F6q,+1%2F24,+-1%2F120%7D,%7B1,+-2%2B4%2F3q,+2%2F3,+-4%2F15%7D,%7B1,+-3%2B9%2F2q,+27%2F8,+-81%2F40%7D%7D%5E(-1)%7B%7Ba%7D,%7Bb%7D,%7Bc%7D,%7Bd%7D%7D

Also with my second result: https://www.wolframalpha.com/input/?i=%7B%7B1,+1,+1%2F2,+1%2F6,+1%2F24,+1%2F120%7D,%7B0,+q,+q,+1%2B1%2F2q,+1%2B1%2F6q,+1%2F2%2B1%2F24q%7D,%7B0,+0,+1,+1,+1%2F2,+1%2F6%7D,%7B1,+0,+0,+0,+0,+0%7D,%7B1,+-1,+1%2F2,+-1%2F6,+1%2F24,+-1%2F120%7D,%7B1,-2,2,-4%2F3,2%2F3,+-4%2F15%7D%7D%5E(-1)%7B%7Ba%7D,%7Bb%7D,%7Bc%7D,%7Bd%7D,%7Be%7D,%7Bf%7D%7D

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')

Answer 2

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}$$