2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You can find a closed form solution that you can use for comparison of the numerical approximations. $\endgroup$ – Amzoti Dec 12 '13 at 23:15
  • $\begingroup$ What is a closed form solution and how would I obtain one? Thanks $\endgroup$ – Levi Dec 12 '13 at 23:25
  • $\begingroup$ Using Undetermined Coefficients, we have $y(x) = (x+1)^2 - e^{x/4 - 1} \csc(\sqrt{7}) \sin(\frac{\sqrt{7}x}{4})$. This is an exact (closed-form) solution. You can use this to compare the numerical results you are getting. $\endgroup$ – Amzoti Dec 13 '13 at 0:52
  • $\begingroup$ Ok I understand. Thank you. In terms of my matrices above, are they correct? $\endgroup$ – Levi Dec 13 '13 at 9:56
2
$\begingroup$

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 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 Finite Difference and Amzoti's solutions picture

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.