1. Heat equation 1D.



2. Wave equation 1D.



3. Optional: Heat and Wave equations 2D and 3D.

1. Introduction.



2. One-dimensional.



3. Multi-dimensional.



4. Finite Element Approximation.

Matlab PDEtoolbox.

Matlab PDEtoolbox documentation.

1. Write your own code to calcuate the 6 FE matrices (M,K,A,Q,G,F) from the mesh quantities P, E, T and the coefficients of the PDE and the boundary condition.
2. Compute the computational time of your code to generate the 6 FE matrices and compare with the computational time of the Matlab code. Use hmax = [0.5,0.2,0.1,0.08]; Maybe hmax = 0.05 if your computer can run a problem of this size.
3. What is the ratio between the time of your code and the time of the Matlab code.
4. How do the computational times and the ratio change as a function of the number of elements?

Specific directions

Show the times to generate your FE matrices and the times Matlab uses to generate the FE matrices.

For timing: use cputime instead of tic toc; run the timing tests 10 times to get the average CPU time. Label the cputime as seconds (cputime returns the CPU time in seconds that has been used by the MATLAB process since MATLAB started). Plot both your time and Matlab time, not just the ratio. On the x-axis, use the number of nodes in the finite elements mesh (size of P). Show that max(max(abs(FEM_yours-FEM_matlab))) is small for all the examples you run.

Note concerning the Matlab 2013 code given below:
the Matlab 2013 code defines the mass matrix as from the term a*u so I set a = 1 to get our mass matrix m_ij = int_Omega phi_i*phi_j (D_COEFF = 1 for us).
To get the "A" matrix, I just did A = A_COEFF*M, but this works only if A_coeff is a constant and does not depend on space variable (x,y).

Matlab 2016 files (Updated Monday Sept 18 at 10:38pm).

Matlab 2013 files (Updated Monday Sept 18 at 11:33pm).

1. Write a Matlab function evaluating the 1 dimensional convolution integral of the initial condition IC(x) and the Green's function on [0,l] for the homogeneous Neumann boundary conditions. Let the user choose the number of terms in the infinite sum, the initial condition, the length of the interval, the evaluation time, and the diffusion coefficient.

2.
A. For the Matlab PDEtoolbox, make the domain a rectangle [0,l]x[0,width], and set the PDE to u_t = (u_xx+u_yy) with homogeneous Neumann boundary conditions. Choose the initial condition (no y dependence) u(x,0) = exp(-(x-x0)^2/w) where w controls the width of the Gaussian and x0 controls of the center of the Gaussian. Compute the approximate solution u_matlab(x,y,t) using the Matlab PDEtoolbox.
B. Compute the convolution integral in 1 dimension to get u_exact(x,y,t) (no dependence in y).
C. Plot ||u_exact(x,y,tfinal)-u_matlab(x,y,tfinal)||/||u_exact(x,y,tfinal)|| in the 2-norm at the FE nodes for several values of hmax. What is the convergence order in h in the 2-norm of the Matlab solution to your exact solution? Suppose error(h) = constant*h^(power), plot log(error) vs log(h) to get the numerically computed order of convergence (the power, h^power).


Specific directions

The domain is a rectangle. The boundaries on the 4 sides of the rectangle. The normal directions to the vertical sides of the rectangle are + and - x direction, the normal to the horizontal sides of the rectangle are the + and - y directions. (In theory u_y, u_yy are zero because your initial condition does not depend on y.)
Remember the nodes of the FE mesh are stored in the variable P. The first row of P contains the x coordinates and the second row of P contains the y coordinates. The number of columns in P is the number of nodes in the FE mesh.

For extra work (only if you have extra time, otherwise, start Project 3), you can do the following: A. Figure out how to determine the number of terms to use in the infinite sum to get a good approximation for a fixed final time tf? B. For a more general code that has dependence in both x and y, you can use the IC(x,y) = alpha(x)*beta(y) and use the product of the one dimensional formulas. C. You can also add the time and space convolution integrals associated with the source term to the PDE. D. Numerically compute {u_x, u_y, u_xx, u_yy, u_t} (using pdegrad or pdecgrad functions for space derivatives and gradient function for the time derivative) for u_exact and u_matlab for several values of hmax and dt in tlist. For each hmax, choose dt small enough so that u_t does not change if you decrease dt. For this dt, plot ||u_t(x,y,tfinal)-(u_xx(x,y,tfinal)+u_yy(x,y,tfinal))|| in the 2-norm for several values of hmax.

1. Use the Matlab PDEtoolbox, solve the wave equation u_tt = (u_xx+u_yy) on a rectangle subject to homogeneous Neumann boundary conditions. Use the initial conditions (no y dependence) u(x,0) = exp(-(x-x0)^2/w) and u_t(x,0) = 0.
2. Use Y = [U; dU/dt] and rewrite the function to pass into the ODE solver so you can solve the wave equation. Compare your PDE solution using the ODE solver with the PDE solution obtained using the Matlab PDEtoolbox.

Specific directions

Change the definition of the mass matrix to pass into the ODE solver.

For extra work (only if you have extra time), you can do the following: Write a Matlab function computing the exact solution to the homogeneous Neumann problem in one dimension using the convolution formulas in Wave equation, convolution formulas. Check the Matlab solution above is close to the exact solution. Use as the domain for the Matlab PDEtoolbox the rectangle.