2019
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.
Preparation for the course
1. Download Matlab and Matlab PDE Toolbox
2. Download SpinDoctor from GitHub
3. Run the examples from SpinDoctor paper
Scientific content of the course
Homework
1. Write a function that computes the stiffness matrix yourself and compare with the Matlab PDE Toolbox stiffness matrix K for c(x,y) = 1.
2. Write a function that gives the forcing term f(x,y) to the PDE.
3. Write a clean program SolvePDE.m that solves
The heat equation d(x,y) du/dt - grad dot (c(x,y) grad u)+a(x,y)u = f(x,y)
subject to
Neumann boundary conditions:
c(x,y) grad u dot n + q(x,y) n = g(x,y)
or Dirichlet boundary conditions:
h(x,y) u = r(x,y);
Allow the geometry to be:
a disk of radius R;
a rectangle of width W and height H.
Write the following as functions depending on User Parameters:
d(x,y)=gg: MassTerm
c(x,y)=gg+cc*exp(-((x-x0)^2+(y-y0)^2)/aa): DiffusionTerm
a(x,y)=gg: AbsorptionTerm
f(x,y)= gg+cc*exp(-((x-x0)^2+(y-y0)^2)/aa): ForcingTerm
q(x,y)=gg: NeumannQTerm
g(x,y)=gg: NeumannGTerm
IC(x,y) = gg+cc*exp(-((x-x0)^2+(y-y0)^2)/aa): InitialCondition
Allow the User to choose
the ode solver (ode45,ode15s,ode23s,ode23t);
the relative and absolute tolerance for the ODE solver;
hmax, the requested finite elements size;
4. In SolvePDE.m, plot the following functions on the finite elements nodes P and plot them using the command pdeplot
DiffusionTerm
ForcingTerm
InitialCondition
Plot each function in its own Figure. Label x and y axes, give a title to the Figure.
5. Solve the ODE system and plot the solution at the time points indicated in UserParameters.m. Plot the solution in a Figure. Label x and y axes, give a title to the Figure. Give the time point of the solution that is displayed. Use 'pause(1)' between the display of successive time points.
6. Make an input file (UserParameters.m) that contains all the necessary user specified parameters for parts 3), 4), 5) above. In your code, if the User chooses a parameter outside of the valid range (such as choosing parms_bc.type = 3) please display a warning and stop the code (disp('Warning, parameter choice not valid. Stopping.'); stop;)
7. Write a function that computes the stiffness matrix yourself and compare with the Matlab PDE Toolbox stiffness matrix K for arbitrary c(x,y).
8. Write a function that computes the G matrix yourself and compare with the Matlab PDE Toolbox matrix G for arbitrary g(x,y).
9. Solve the wave equation.
10. Integrate the heat PDE over Omega and write down the expression for d/dt (int_Omega u(x,y,t) dxdy), using Green's identities. Numerically compute int_Omega u(x,y,t) dxdy and compare with int_s int_Omega f(x,y,s) dxdy dt.
The most recent version of UserParameters.m
Turn in the homework by making a folder called HW1, HW2, etc and put all necessary files inside each folder. Send the folder direct to me by email.
HW1. Write a function that computes the stiffness matrix yourself and compare with the Matlab PDE Toolbox stiffness matrix K for c(x,y) = 1.
HW2. Write a function that gives the forcing term f(x,y) to the PDE.
HW3. Do 3), 4), 5), 6) above.
Turn in UserParameters.m as the file contains all the user input parameters.
Turn in SolvePDE.m to solve the heat equation (and any other needed files to run the code).
At the end of SolvePDE.m put your function definitions for PDE (give the functions the indicated names):
MassTerm
DiffusionTerm
AbsorptionTerm
ForcingTerm
NeumannQTerm
NeumannGTerm
InitialCondition
I will be changing the UserParameters.m file to test that your code works in all choices of User Parameters.
HW4. Write a function that computes the stiffness matrix yourself and compare with the Matlab PDE Toolbox stiffness matrix K for arbitrary c(x,y).
I will change the function DiffusionTerm to test your code.
HW5. Write a function that computes the G matrix yourself and compare with the Matlab PDE Toolbox matrix G for arbitrary g(x,y).
I will change the function NeumannGTerm to test your code.
HW6. Solve the wave equation. Turn in UserParameters.m as the file contains all the user input parameters.
Turn in SolvePDE.m to solve the heat equation (and any other needed files to run the code).
HW7. Integrate the heat PDE over Omega and write down the expression for d/dt (int_Omega u(x,y,t) dxdy), using Green's identities. Numerically compute int_Omega u(x,y,t) dxdy and compare with int_s int_Omega f(x,y,s) dxdy dt.
Software used for the class
Optional additional reading
Jing-Rebecca Li
Last modified: Thu Oct 26 16:22:44 GMT 2017