Trouble in using finite difference method to solve a boundary value problem (attempts and pictures included)

finite differencesfinite element methodmatricesnumerical methodsordinary differential equations

I want to numerically solve (using FDM)$$-y''(t)+2y'(t)=1, t\in (0,1)\\y(0)=1, y(1)=3$$

First, I check with symbolab that the analytic solution would be $1-\frac{3}{2\left(-1+e^2\right)}+\frac{3}{2\left(-1+e^2\right)}e^{2t}+\frac{t}{2}$:

$\frac{d}{dt}\left(1-\frac{3}{2\left(-1+e^2\right)}+\frac{3}{2\left(-1+e^2\right)}e^{2t}+\frac{t}{2}\right) = \frac{3}{e^2-1}e^{2t}+\frac{1}{2}$, and

$\frac{d^2}{dt^2}\left(1-\frac{3}{2\left(-1+e^2\right)}+\frac{3}{2\left(-1+e^2\right)}e^{2t}+\frac{t}{2}\right) = \frac{6}{e^2-1}e^{2t}$,

and we can easily check that the boundary conditions are satisfied.

Now, I attempt this problem numerically:
Using the centered difference approximations for the first and second derivative, I get two equations:

$D^2y_j = (y_{j-1}-2y_j+y_{j+1})/h^2\\Dy_j = (y_{j-1}+y_{j+1})/2h$.

I use the second order equation of the problem to derive the relation $(\frac{-1}{h^2} + \frac{1}{h})y_{j-1}+\frac{2}{h^2}y_j+(\frac{-1}{h^2} + \frac{1}{h})y_{j+1}$.

For example, if we take a mesh size of $.25$, we would have to solve 3 unknowns at $.25, .5, .75$, and our matrix equation would look like $$\begin{pmatrix}32&-12&0\\ -12&32&-12\\ 0&-12&32\end{pmatrix}v = \begin{pmatrix}1+12\\ \:\:1\\ \:\:1+12\cdot 3\end{pmatrix} = \begin{pmatrix}13\\ 1\\ 37\end{pmatrix}\\v= \begin{pmatrix}y_{.25}\\ \:\:y_{.5}\\ y_{.75}\end{pmatrix}$$

But if we solve this $\begin{pmatrix}y_{.25}\\ \:\:y_{.5}\\ y_{.75}\end{pmatrix} = \begin{pmatrix}\frac{67}{92}\\ \frac{79}{92}\\ \frac{34}{23}\end{pmatrix}$, and plot it, it looks off from the analytic solution.

One would think that increasing the number of nodes between 0 and 1 would be better, but increasing the unknown to 500 only made it worse:

Now, I may have made the program incorrectly, but the matrix vector equations seem to come out fine. I will upload the code if anyone is interested, but I think that my problem comes from a misunderstanding of the algorithm of the finite difference method. Any help would be lovely.

Best Answer

$$\color{gray}{Dy_j=\frac{y_{j−1}+y_{j+1}}{2h}}$$ is wrong and is in this wrong form used in the further equations. The central difference quotient is $$ Dy_j=\frac{y_{j+1}-y_{j-1}}{2h}=\frac{-y_{j−1}+y_{j+1}}{2h}. $$ It might be that the minus sign in the last form was simply lost sometime.

Related Solutions

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

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

[Math] Finite Difference Boundary Conditions

With Dirichlet conditions you can either build it into the right hand side of the equation, or you can introduce a trivial equation for those conditions as actual variables. It will be the same.

With Neumann or Neumann-like (e.g. Robin) conditions you want to be careful that you discretize the boundary derivative in a way which has the same order as your discretization of the interior derivatives. One way to do this with finite differences is to use "ghost points". I confess that this is rather hard to motivate within the finite difference framework but it gives results that are much like those you get in the finite element framework.

The point here is that the values $u_0,u_{n+1}$ are thought of as being variables, and there is an equation involving the "interior" derivative at these points. You then add additional points $u_{-1},u_{n+2}$ and choose their values to satisfy the derivative condition. For example, with homogeneous Neumann conditions you take the approximation of the derivative at the first point to be $\frac{u_1-u_{-1}}{2\Delta}$ and set that to zero so that $u_{-1}=u_1$. Same for $u_{n+2}$. This then actually enters into the problem because you have an equation for the second derivative at the leftmost point: $u_{-1}-2u_0+u_1=b_0$, and $u_{-1}=u_1$ so this is $-2u_0+2u_1=b_0$.

Best Answer

Related Solutions

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

[Math] Finite Difference Boundary Conditions

Related Question