Linear systems that arise in many applications can become quite large. It is often necessary to exploit any structure and /or sparsity in the matrices to reduce the computation burden. We will consider the one-dimensional Poisson equation: -y"(t) = f(t) on [0, 1], with y(0) = y(1) = 0.

Such problems are frequently very hard to handle because it is often not possible to express y(t) in terms of elementary functions (as is done in undergraduate courses on ODE). Numerical methods are employed in order to approximate y(t) at discrete points inside the interval [a, b]. This approach will leads to a linear system of equations.

The approach we consider begins by subdividing the interval [a, b] into n + 1 equal subintervals, each of length b-a/n+1: t₀ = a, t₁ = a + h, t₂ = a + 2h, ··· , t_n = a + nh, t_n+1 = b; The points t_i = a + ih are called grid points, and the value h = b-a/n+1 is called the step size. Smaller step sizes generally produce better approximations to the derivatives, so better accuracy requires smaller step size, and hence larger number of grid points.

Let y_i = y(t_i) and f_i = f(t_i), and show that by using the finite difference method, we can compute approximations to y_i by solving the linear system T_y = h²f, where

(a) Prove that the matrix T has an LU factorization with L(i, i) = 1 and L(I + 1, i) = -i/i+1, U(i, i) = i+1/i and U(i, i + 1) = -1.

(b) Is partial pivoting needed?

(c) Write a Matlab function T = poissonmat(n) which the n x n matrix T for given n.

(d) Write a Matlab function y = poissonsolve(f) which, given a vector f of length n, computes the solution y of T_y = h²f.

(e) Let f(t) = sin πt, then the right hand side f_i = sin π(t_i). Experiment with various values of n, say n = 10, 100, 500, 1000. Fill the following table: