Converting a PDE into a matrix equation of the form Ax=b and discuss error estimate using discretization method

354 Views Asked by At

Let $$-u_{xx}+u=f(x)=(\pi-x)\sin x, \ x\in (0,\pi),\ u(0)=0,\ u(\pi)=2.$$ Solve the problem on $$h=\frac{\pi}{N}, N=16,32, 64, 128, 256$$ $$0=x_{o}<x_{1} < \dots <X_{N}=\pi,$$ $$x_{I+1}=x_{1}+h, I=0,\dots N-1. $$

Question 1: Convert the problem to $Au=b$. $A=?$, $b=?$.

Question 2: Discuss about the error estimate

Question 3: Solve $Au=b$ by a direct solver. Plots. At $x=\frac{\pi}{2}$, repeat the solution of each grid.

What I know is that I want to use the the Taylor expansion of $f(x+h)$ up to order $2$. This I know the expression but then I am not sure how to use this information to find matrix $A$ and vector $b$. This is numerical analysis and I am not sure how to go about answering these questions. Any help with some explanations would really be much appreciated.

1

There are 1 best solutions below

1
On

This is a linear ODE, which we can rewrite as $$u''(x) = (x-\pi)\sin(x) + u(x)$$

Then, we have that $$\frac{1}{h^2}(u(x_{i-1})-2u(x_i)+u(x_{i+1})) = (x_i-\pi)\sin(x_i) + u(x_i)$$ for $i = 2,\dots , n-2$, using a central difference to approximate $u''(x_i)$. For $i=1$, we have $u(x_{i-1}) = u(0) = 0$ and for $i = n-1$, we have $u(x_{n-1+1}) = u(\pi) = 2$.

Using this, we rearrange the equations:

For $i=1$, $$ (2+h^2)u(x_1)-u(x_2) = h^2(\pi-x_1)\sin(x_1) $$

For $i = 2,\dots, n-2$, $$ -u(x_{i-1})+(2+h^2)u(x_i)-u(x_{i+1}) = h^2(\pi-x_i)\sin(x_i) $$

For $i = n-1$, $$ -u(x_{n-2})+(2+h^2)u(x_{n-1}) = h^2(\pi-x_{n-1})\sin(x_{n-1}) + 2 $$

This is a system of linear equations, which can be written as a (tridiagonal) matrix vector equation. The right hand sides constitute the $\mathbf b$ vector, while the coefficients of the left hand sides define the $A$ matrix. Of course, $u_i = u(x_i)$ defines the $\mathbf u$ vector. We can write the system as $A\mathbf u = \mathbf b$.


We can compute an error bound for linear systems $A\mathbf u=\mathbf b$ by the following inequality

$$ E(\mathbf u) = \frac{\|\hat{\mathbf u} - \mathbf u\|}{\|\mathbf u\|} \leq \text{cond}(A)\frac{\|\hat{\mathbf b} - \mathbf b\|}{\|\mathbf b\|} = \text{cond}(A)E(\mathbf b)$$

where $\text{cond}(A) = \|A\|\|A^{-1}\|$ represents the "closeness" of $A$ to a singular matrix. The condition number can be estimated with the upper bound $\text{cond}(A) \leq 1+4h^{-2}$, calculated using this method. Thus, the (rough) error bound for computing $\mathbf u$ is given by $$ E(\mathbf u) \leq (1+4h^{-2})E(\mathbf b) $$