Scientific Computing, homework assigned Mar 10/12, due in class Mar 24. Part 2. Modify the code of homework 5 (Van Loan, exercise 7.3.1) to use the alternative discretization discussed in class. There is a snag: the modified discretization is not symmetric, so one must use tridiagonal L*U decomposition instead of tridiagonal Cholesky. Van Loan supplies codes, and I've included them in the Matlab context (see the link further below). The modified discretization is obtained as follows. For each interior point xi, consider the "tent" function, tenti(x), which is defined as the piecewise linear function that is 0 for x<=xi-1, 1 for x=xi, and 0 for x>=xi+1. It increases linearly from 0 to 1 on the interval (xi-1,xi) and decreases linearly from 1 to 0 on the interval (xi,xi+1). Take the differential equation, multiply it by tenti, and integrate from xi-1 to xi+1. So: ( Int x : xi-1<=x<=xi+1 : -D[p(x)Du(x)]*tenti(x)*dx ) + ( Int x : xi-1<=x<=xi+1 : q(x)*u(x)*tenti(x)*dx ) = ( Int x : xi-1<=x<=xi+1 : r(x)*tenti(x)*dx ) (Notation: (Int x : range(x) : f(x)*dx) means the integral of f with respect to x over the indicated range.) For the first term we use integration by parts, noting that tenti vanishes at the endpoints. So ( Int x : xi-1<=x<=xi+1 : -D[p(x)Du(x)]*tenti(x)*dx ) = ( Int x : xi-1<=x<=xi+1 : p(x)*Du(x)*Dtenti(x)*dx ) Dtenti(x) is piecewise constant; it has the value 1/h on (xi-1,xi) and the value -1/h on (xi,xi+1). So ( Int x : xi-1<=x<=xi+1 : -D[p(x)Du(x)]*tenti(x)*dx ) = ( Int x : xi-1<=x<=xi+1 : p(x)*Du(x)*Dtenti(x)*dx ) = (1/h)*(Int x: xi-1<=x<=xi: p(x)*Du(x)*dx) - (1/h)*(Int x: xi<=x<=xi+1: p(x)*Du(x)*dx) Now, if p(x) is a constant then the result of these integrals just involves u at the meshpoints. In the general case we need to approximate the integrals in terms of u at the mesh points. In class I suggested the approximation: (Int x: xi<=x<=xi+1: p(x)*Du(x)*dx) = (1/2)*(p(xi)+p(xi+1))*(u(xi+1)-u(xi)) but another possibility (which turns out to be somewhat more accurate) is (Int x: xi<=x<=xi+1: p(x)*Du(x)*dx) = p(xi+h/2)*(u(xi+1)-u(xi)) This approximation also makes the whole process closer to the original exercise 7.3.1, so I suggest that you use it. In addition we need approximations for the other two terms: ( Int x : xi-1<=x<=xi+1 : q(x)*u(x)*tenti(x)*dx ) and ( Int x : xi-1<=x<=xi+1 : r(x)*tenti(x)*dx ) To obtain these approximations one does the standard thing: construct weights for the points xi-1, xi and xi+1 so that a three-point approximation is exact when (q(x)*u(x)) or r(x) is a polynomial of at most 2nd degree. The correct weights turn out to be h*(1,10,1)/12, so we obtain the approximations: ( Int x : xi-1<=x<=xi+1 : q(x)*u(x)*tenti(x)*dx ) = (h/12)*(q(xi-1)*u(xi-1)+10*q(xi)*u(xi)+q(xi+1)*u(xi+1)) and ( Int x : xi-1<=x<=xi+1 : r(x)*tenti(x)*dx ) = (h/12)*(r(xi-1)+10*r(xi)+r(xi+1)) With this you have all the expressions needed to write down the three-point discretization of the equation around the interior point xi. Observe that the discretization is not exactly symmetric (it would be symmetric if q(x) is a constant). So to solve it we should use tridiagonal L*U decomposition. Van Loan discusses this in Section 6.2.1 and supplies codes to carry out the decomposition and the solution of the resulting pair of bi-diagonal equations. I've prepared all this in the Matlab script that follows below. Matlab context. This script contains almost a complete running code for the exercise, and it is a modification of my posted solution to homework 5. The only thing that is left out of the script is the section that discretizes the interior equations; that section is marked by a sequence of question marks. Your task is to fill in the missing section and then run the script. Please hand in your code and output. If you do it correctly I expect the following output: >> twopoint1 original function p of 7.3.1 n u(0.25) u(0.50) u(0.75) du(0.50) 5 -3.6545781e-01 -3.2439023e-02 2.9284310e-01 -3.24e-02 9 -4.2535487e-01 -2.0092434e-02 3.8951266e-01 1.23e-02 17 -4.6423349e-01 -1.0166879e-01 2.8919203e-01 -8.16e-02 33 -3.6327434e-01 -1.0504016e-02 3.3579685e-01 9.12e-02 65 -3.8459630e-01 -3.1062642e-02 3.2161378e-01 -2.06e-02 129 -3.8531546e-01 -3.1625220e-02 3.2126606e-01 -5.63e-04 257 -3.8532940e-01 -3.1655214e-02 3.2122684e-01 -3.00e-05 513 -3.8533296e-01 -3.1662767e-02 3.2121708e-01 -7.55e-06 1025 -3.8533386e-01 -3.1664658e-02 3.2121464e-01 -1.89e-06 modified function pmod, constant n u(0.25) u(0.50) u(0.75) du(0.50) 5 -3.6545781e-01 -3.2439023e-02 2.9284310e-01 -3.24e-02 9 -4.4207284e-01 -8.8217801e-02 2.7478046e-01 -5.58e-02 17 -4.6420565e-01 -1.0715264e-01 2.6465621e-01 -1.89e-02 33 -3.6687740e-01 -2.3936962e-02 3.0397896e-01 8.32e-02 65 -3.8540689e-01 -3.8870187e-02 2.9693888e-01 -1.49e-02 129 -3.8603891e-01 -3.9294773e-02 2.9673882e-01 -4.25e-04 257 -3.8603838e-01 -3.9294159e-02 2.9673912e-01 6.14e-07 513 -3.8603836e-01 -3.9294126e-02 2.9673914e-01 3.26e-08 1025 -3.8603836e-01 -3.9294124e-02 2.9673914e-01 1.96e-09 The last column shows the change in the computed value of the solution at the point x=0.50 upon each doubling of the mesh. Pure 2nd order convergence would be shown by having that change go down by a factor of 4 with each mesh doubling, and 4th order convergence would have it go down by a factor of 16.