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.