Homework 5

Due to the dropbox by 11:00pm on November 22, 2013.

Problem 1

The bottom figure on page 172 of ATAP shows the approximate solution to the Airy equation (21.4) computed using a Chebyshev spectral method. The resulting Chebfun has length \(n=184\).

Write a program that computes the approximate solution for values of n = 10:10:200 and produces a semilog plot of the error vs. \(n\). To measure the error, take the vector \(u_n\) obtained with \(n\) points and use Chebfun to generate the interpolating polynomial \(p_n(x)\) over the interval. Also define \(\hat p(x)\) to be the Chebfun approximation of the true solution (allowing it to use as many points as needed to get machine precision). Then define the error as \(E_n = \|p_n-\hat p\|_\infty\) (which can be computed with norm(pn-phat) if these are defined properly as Chebfuns).

Also plot the computed and true solutions together on the same plot for selected values of \(n = 100, 110, 120\).

Note the following:

• The matlab function airy can be used, but you will have to rescale the value returned by this function so that \(\hat p(x)\) at the left end point has the value 1, since this is the boundary condition imposed on the numerical solution.
• Set up and solve an \(n\times n\) linear system for each value of \(n\). You can get the matrix by evaluating L(n) where L is the chebop defined on page 172, and then you will have to define the right hand side of the system appropriately to impose the desired boundary conditions

Problem 2

To see how convergence of the spectral method compares to convergence of the interpolation problem of approximating the Airy function, write a program that uses Chebfun to compute the polynomial \(p_n(x)\) that interpolates at \(n\) Chebyshev points and plot the error \(E_n = \|p_n-\hat p\|_\infty\) vs. \(n\) on the same semilog plot as in figure 1. Again do this on the interval \([-30,30]\) and scale the function you are interpolating to have \(\hat p(-30) = 1\) so that \(\hat p\) is the same as before.

Problem 3

In lecture on Monday Nov. 18 we will investigate a Chebyshev spectral method to solve the advection equation \(u_t + u_x = 0\) using the Forward Euler method for time stepping. This code can be found in the “Spectral-Advection” notebooks listed in Sample codes.

Modify the code to instead use the 4th order Runge-Kutta method for time stepping and investigate what time step restriction must be imposed to obtain a stable approximation, by plotting the eigenvalues of the spectral differentiation matrix being used and taking into account the stability region of the method, which can be plotted by adapting the code in Program 25 of SMM. Confirm this by running your code for various values of \(\Delta t\) to see when it is stable.

The Runge-Kutta method for \(u'(t) = f(u)\) is:

\(k_1 = f(U^n)\)

\(k_2 = f(U^n + 0.5 \Delta t k_1)\)

\(k_3 = f(U^n + 0.5 \Delta t k_2)\)

\(k_4 = f(U^n + \Delta t k_3)\)

\(U^{n+1} = U^n + \frac{\Delta t}{6} (k_1 + 2k_2 + 2k_3 + k_4)\)

Explore how the error at time \(t=1.0\) behaves both as the number of grid points \(n\) is varied and as the time step \(\Delta t\) is varied.

Problem 4

Develop a Chebyshev spectral method to solve the heat equation \(u_t = u_{xx}\) on the interval \([-1,1]\) with boundary conditions \(u(-1)=u(1)=0\). Use Chebfun to determine the spatial differentiation matrix and the implicit Trapezoid rule for time stepping.

Apply this method with initial data \(u(x,0) = \exp(-50 x^2)\) and for which the solution is essentially a decaying Gaussian (the boundary conditions aren’t quite right but for small \(t\) the values of the Gaussian below the rounding error level at the boundaries).

Explore how the error at time \(t=0.02\) behaves both as the number of grid points \(n\) is varied and as the time step \(\Delta t\) is varied.

For a solution to this problem: see the Spectral-Heat notebook on the Sample codes page.