Class page for "Mathematical and numerical foundations of modeling and simulation using partial differential equations"
Part of the Program: French-Vietnam Master in Applied Mathematics

I will update this page regularly until the start of the class. My lectures will be simpler and shorter than the lecture notes provided below.

Background material (I will do a brief review at the beginning of the class).
Scientific content of the class.

Software that is useful to have on your laptop if possible (otherwise I can provide some licenses you can share)
Lectures and projects.
For each project, please create a README file giving directions on how to use your code. What commands to type, etc.

Inside your code, please use "disp('information for the user')" to tell the user what your code is doing at the important steps. Display the results and conclusions as the user runs the code. Put a "pause" between the important steps to let the user "press return" when she or he is ready to see the next results. For example, "disp('Now generating the FE matrices')", "disp('Press return to see the timing tests')".

Use both Figures and the display text in the Matlab window to explain what you are doing for each project.

Each Figure should have a title, xlabel, ylabel, legends (if you have several curves on the same plot).

Lecture 1.

Project 1

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).

Lecture 2.

Project 2

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.

Wave equation, general properties.
Wave equation, convolution formulas.

Project 3

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.

Grade out of 20 points.

Please come to class On Friday Sept 29 at 8AM

On Friday Sept 29 at 11AM

A. Turn in your code for Projects 1, 2, 3 to me. I will check for correctness of and clarity of code. (Count for 8 points)

On Friday Sept 29 at 11AM

B. Turn in the following Matlab Figures with the requested information. (Count for 6 points)

Project 1.

Run your code on the domain [0,1]x[0,0.4] for at least 4 different values of hmax.

Fig 1.Plot the cputime of Matlab and the cputime of your way of generating the FE matrices as a function of the number of elements, N_elements, for the different meshes. The x-axis should be log10(N_elements) for the different meshes. The y-axis should be log10(cputime_matlab) and log10(cputime_yourcode). The legend should say 'cputime matlab' and 'cputime my code'. Do a linear fit of log10(N_elements) vs log10(cputime) and show fitting equation. On the Figure window, click "Tools", "Basic Fitting", check "linear" and "show equations". In the title of the Figure, describe how the cputime behaves as a function of N_elements for Matlab and for your way.

Example Figure 1: the cputime of Matlab PDE toolbox solver.

Fig 2. For the 4 or more different meshes above, plot the difference between your FE matrices and the Matlab FE matrices. Show the differences of all the FE matrices, M,K,A,G,Q,F.

Project 2.

On 4 different FE meshes (4 different values of hmax) solve dt = u_xx+u_yy on the domain [0,1]x[0,0.4] subject to initial condition u(x,y,0) = exp(-(x-0.3)^2/0.01) (no dependence on y) for t in [0,tfinal=0.05] and homogeneous Neumann boundary conditions, using the Matlab PDEtoolbox. Call the solution u_matlab(x,y,t).

Compute the convolution integral of the initial condition with the Green's function to get u_exact(x,y,t).

Fig 3. For the finest mesh, show plot of u_matlab(x,y,0), u_matlab(x,y,tfinal) and u_exact(x,y,tfinal); Use the command
pdeplot(mymodel,'XYData',u,'Contour','on','ColorMap','jet');

Example Figure 3.

Fig 4. Let err = ||u_matlab(:,:,tfinal)-u_exact(:,:,tfinal)||/||u_exact(:,:,tfinal)||. Plot log10(hmax_vec) vs log10(err). Do a linear fit of log10(hmax_vec) vs log10(err) and show fitting equation. On the Figure window, click "Tools", "Basic Fitting", check "linear" and "show equations". Choose your hmax so the 4 data points on this figure are fit well by the linear fit. In the title, give the order of convergence between u_matlab and u_exact in hmax.

Example Figure 4.

Add the following changes in the tolerances of the Matlab PDEtoolbox to your code

%%% Set the tolerances for the solver options
model.SolverOptions.AbsoluteTolerance=1e-7;
model.SolverOptions.RelativeTolerance=1e-5;
model.SolverOptions.ResidualTolerance=1e-5;
mymodel.SolverOptions.ReportStatistics='on';
% solve PDE
R = solvepde(model,tlist);

Fig 5. Replot log10(hmax_vec) vs log10(err). Do a linear fit of log10(hmax_vec) vs log10(err) and show fitting equation. In the title, give the order of convergence between u_matlab and u_exact in hmax.

Example Figure 5.

Project 3.

On 4 different FE meshes (4 different values of hmax) solve d_tt = u_xx+u_yy on the domain [0,1]x[0,0.4] subject to initial condition u(x,y,0) = exp(-(x-0.3)^2/0.01) (no dependence on y) and u_t(x,y,0) = 0 for t in [0,tfinal=0.5] and homogeneous Neumann boundary conditions, using the Matlab PDEtoolbox. Call the solution u_matlab(x,y,t).

Let utilde(:,t) = [u(:,t); u_t(:,t)], use the Matlab ODE solver on utilde(:,t) with a modified mass matrix and a modified right hand side function to obtain u_odesolver(x,y,t).

Fig 6. Plot u_matlab(x,y,0), u_matlab(x,y,tfinal) and u_odesolver(x,y,tfinal). Use the command
pdeplot(mymodel,'XYData',u,'Contour','on','ColorMap','jet');

Example Figure 6.

On Friday Sept 29 starting 12pm

Oral exam, 10 minutes per person. (Count for 6 points)

Extra work (Count for maximum 2 points)

Jing-Rebecca Li
Last modified: Thu Oct 26 16:22:44 GMT 2017