[Math] Solve the following boundary value problem using the finite difference method.

Question

Solve $$y''=\frac{1}{2}y'-\frac{1}{2}y+\frac{x^2+3}{2}, ~~~~~y(0)=1, ~~y(4)=24$$
using the second order finite difference approximation order with $h=1$.

I know that we use $y''=p(x)y'(x)+q(x)y(x)+r(x)$, and here $p(x)=\frac{1}{2}$, $q(x)=-\frac{1}{2}$ and $r(x)=\frac{x^2+3}{2}$, and to solve this problem we need to create matrices of the form $Aw=b$.

Our $A$ matrix is of the form:
$$
\begin{bmatrix}
2+h^2\left(\frac{1}{2}(x_1)\right) & -\left(1-\frac{h}{2}(-\frac{1}{2}(x_1))\right) & 0 & & 0 & 0 \\
-\left(1+\frac{h}{2}(-\frac{1}{2}(x_2))\right) & 2+h^2\left(\frac{1}{2}(x_2)\right) & -\left(1-\frac{h}{2}(-\frac{1}{2}(x_1))\right)& & & 0 \\
0 & \ddots & \ddots & \ddots & \ddots & \\
& \ddots & \ddots & \ddots & \ddots & 0 \\
0 & & \ddots & \ddots & \ddots & -\left(1-\frac{h}{2}(-\frac{1}{2}(x_1))\right) \\
0 & 0 & & 0 & -\left(1+\frac{h}{2}(-\frac{1}{2}(x_n))\right) & 2+h^2\left(\frac{1}{2}(x_n)\right)
\end{bmatrix}
$$

and this simplifies to
$$\begin{bmatrix}
2.5 & -1.25 & 0 & & 0 & 0 \\
1.25 & 2.5 & -1.25 & \ddots & & 0 \\
0 & \ddots & \ddots & \ddots & \ddots & \\
& \ddots & \ddots & \ddots & \ddots & 0 \\
0 & & \ddots & \ddots & \ddots & -1.25 \\
0 & 0 & & 0 & 1.25 & 2.5
\end{bmatrix}.$$

Our $b$ is
$$\begin{pmatrix}
-\frac{x_1^2+3}{2}+\frac{3}{4}w_0 \\
-\frac{x_2^2+3}{2} \\
\vdots \\
-\frac{x_{n-1}^2+3}{2} \\
-\frac{x_n^2+3}{2}+1.25w_{n+1}
\end{pmatrix}.$$

Then we solve $Aw=b$ by using LU-decomposition.

However, I am very doubtful that my matrices are correct and I am really struggling with finding a useful and clear example that illustrates this problem. If anyone can help me solve this I would be grateful. Also, if anyone knows how to code this on MATLAB I would appreciate it.

It would be interesting to see the plots to this.

Thank you.

Answer 1

This is a MATLAB realization of your problem with matrix you suggest.

    %% ;
k = 5; % number of nodes of the grid
D = sparse(1:k,1:k,2*ones(1,k),k,k); %diagonal
E = sparse(2:k,1:k-1,-1.25*ones(1,k-1),k,k); %bottom diagonal
A = (-1)*E + D + E'; %E'  transposed bottom diagonal, i.e. upper one
x = linspace(0,4,k); % vector of independent variable
w_bound = zeros(1,k); w_bound(1) = 3/4; w_bound(k) = 1.25 * 24; % it's your w0 and wn I don't understant how you get it
b = -0.5 * (x'.^2 + 3) + w_bound'; % the right part of matrix equation
%% Constructing picture 
figure(1)
hold on
plot(x, A \ b);
x_a = linspace(0,4,100);
y_a = (x_a + 1).^2 - exp(0.25 * x_a - 1) .* (1 / sin(sqrt(7))) .* sin(0.25 * sqrt(7) * x_a);
plot(x_a, y_a);
legend('your soulution', 'Amzoti''s Solution')

Picture with results of comparison your and Amzoti's solutions

It seems that your matrix A is not correct. Let's consider it again more accurately.

We have $$y''=\frac{1}{2}y'-\frac{1}{2}y+\frac{x^2+3}{2}, ~~~~~y(0)=1, ~~y(4)=24,$$ and $$y''=\frac{y_{i-1}-2y_i+y_{i+1}}{h^2}, y' = \frac{y_{i+1}-y_{i-1}}{2h}, y_0 = 1, y_n =24.$$ I used here central difference for the first derivative. Let's sum up it together and collect terms:

$$ \frac{y_{i-1}-2y_i+y_{i+1}}{h^2} - \frac{1}{2}\frac{y_{i+1}-y_{i-1}}{2h} + \frac{1}{2}y_i = \frac{x_i^2+3}{2},y_0 = 1,y_n = 24$$, or

$$y_{i-1}\left(\frac{1}{h^2} + \frac{1}{4h}\right) + y_i\left(\frac{-2}{h^2}+\frac{1}{2}\right) + y_{i+1}\left(\frac{1}{h^2}-\frac{1}{4h}\right) = \frac{x_i^2+3}{2},y_0 = 1,y_n = 24$$.

I used the same $h = 1$ as you did.

MATLAB code:

%% the same steps only for my procedure
k       = 5;
x       = linspace(0, 4, k);
h       = (x(end) - x(1)) / (k - 1);
D       = sparse(1 : k, 1 : k,(-2 / h^2 + 0.5) * ones(1,k), k, k);
D_low   = sparse(2 : k, 1 : k-1,(1  / h^2 + 1 / 4 /h) * ones(1, k - 1), k, k);
D_up    = sparse(1 : k - 1, 2 : k,(1  / h^2 - 1 / 4 / h) * ones(1, k - 1), k, k);
A       = D + D_low + D_up;
A(1,1) = 1; A(1,2) = 0;
A(end,end) = 1; A(end,end - 1) = 0;
b = (0.5 * (x.^2 + 3))'; b(1) = 1; b(end) = 24;

%% Pictures for finite difference solution
figure(2)
hold on
plot(x, A^(-1) * b);
x_a = linspace(0,4,100);
y_a = (x_a + 1).^2 - exp(0.25 * x_a - 1) .* (1 / sin(sqrt(7))) .* sin(0.25 * sqrt(7) * x_a);
plot(x_a, y_a);
legend('Finite difference method','Amzoti''s solution');

The matrix $A$ is: $$\begin{bmatrix} 1 & 0 & 0 & & 0 & 0 \\ 1.25 & -1.5 & 0.75 & \ddots & & 0 \\ 0 & \ddots & \ddots & \ddots & \ddots & \\ & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & 1.25 & -1.5 & 0.75 \\ 0 & 0 & & 0 & 0 & 1 \end{bmatrix}.$$ Finite Difference and Amzoti's solutions picture