CMPSC/Mathematics 451
MATLAB Program One
Due 6 October 2021
Fall 2021
This MATLAB assignment has two parts.
Part One. Write a MATLAB/Octave function that solves the linear
system
Av = gvec (1)
with a symmetric, positive definite, tridiagonal m×m matrix A given by
A =
a1 b1 · · · · · · · · · 0
b1 a2 b2 · · · · · · · · ·
. . .
. . .
. . .
. . .
. . .
. . .
bm−2 am−1 bm−1
bm−1 bm
.
The diagonal entries are stored in a = (a1, . . . , am)
T and the off-diagonal
entries are stored in b = (b1, . . . , bm−1)
T . The first line of the function will
be
function v =thomas chol(a,b,g vec)
which delivers the solution to (4).
Test your function with the MATLAB lines
� b=randn(10,1);
� a=2*abs(b);
� a=[a; max(a)];
� T=diag(a)+diag(b,-1)+diag(b,1);
� g vec=T*ones(11,1);
� v=thomas chol(a,b,g vec);
Your answer for v should be all ones or very close to it. If you wish, you
can also compare your result to that from vhat=T\g_vec.
1
Part Two. Write a MATLAB/Octave function to solve the two-point
boundary value problem
−v′′(t) = g(t), 0 < t < 1 (2) v(0) = vleft, v(1) = vright (3) This problem is based upon Example 4.17 on p.85-87 of Ascher and Grief. As explained in class, you construct the linear system Av̂ = gvec (4) where gvec = h2g(t1) + vleft h2g(t2) ... h2g(tN−2) h2g(tN−1) + vright h = 1/N, tj = jh (5) and A is the symmetric, positive definite tridiagonal matrix A = 2 −1 · · · · · · · · · 0 −1 2 −1 · · · · · · · · · . . . . . . . . . . . . . . . . . . −1 2 −1 −1 2 . (6) The resulting vector v = vleftv̂ vright . yields the approximate solution v(tj) ≈ vj, j = 0, . . . , N . This second function will have the first line function v =diffeq setup(g, v left, v right,N) where g is a function in (2), v left and v right correspond to the left and right boundary condition in (3) and N is the number of mesh points. The function thomas chol performs the following tasks 2 1. Compute h = 1/N , set up mesh t = (0:h: 1); 2. Compute g vec above from the function g,v left, and v right as in (5) 3. Set a = 2 ∗ ones(N − 1, 1) and b = −ones(N − 2, 1) to implicitly construct A in (6) 4. Solve (4) using thomas solve to obtain vector vhat. 5. Set v = [v left; vhat; v right]; 6. Plot v against t with appropriate labels. I am going to provide you with three functions in the MATLAB files g1.m, g2.m and g3.m on Canvas. For the functions g1 and g2 the boundary conditions will be vleft = vright = 0. For g3, the boundary conditions will be vleft = 0, vright = 16/( √ 2 ∗ π2). In all cases, set N = 100. I strongly rec- ommend writing a short MATLAB/Octave script which calls diffeq setup three times with the appropriate parameters. Octave Note. When using a function as a argument in Octave, it is best to call it with the line v = diffeq setup(@g, v left, v left,N) The @ is necessary. 3