Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This tour explores the numerical computation of PDE on meshes. The spacial derivatives are replaced by discrete Laplacian operators on meshes. The time derivatives are computed with explicit time stepping.
from __future__ import division
import nt_toolbox as nt
from nt_solutions import meshproc_5_pde as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
We first define differential operator using sparse matrices.
Load a mesh of $n$ vertices.
name = 'bunny'
clear options
options.name = name
[X0, F] = read_mesh(name)
n = size(X0, 2)
options.name = name
Display it.
plot_mesh(X0, F, options)
Exercise 1
Compute the cotangent Laplacian $W$ of the mesh $$ W_{i,j} = \text{cot}(\al_{i,j}) + \text{cot}(\be_{i,j}) $$ where $\al_{i,j}$ and $\be_{i,j}$ are the two angles oposite the the edge $(i,j)$.
solutions.exo1()
## Insert your code here.
Compute the symmetric Laplacian matrix $L=D-W$ where $D = \text{diag}_i( \sum_j W_{i,j} )$
d = full(sum(W, 1))
D = spdiags(d(: ), 0, n, n)
L = D - W
Compute normalized operators $\tilde W = D^{-1}W$ and $\tilde L = D^{-1}L$.
iD = spdiags(d(: ).^(-1), 0, n, n)
tW = iD * W
tL = iD * L
The heat equation is the prototype of parabolic PDE. It produces a smoothing of the initial data that increases with time.
This PDE computes a function $f(x,t)$ of space and time that solves $$ \pd{f}{t} = -\tilde L f \qandq f(\cdot,t=0) = f_0. $$
Load a texture image.
M = load_image('lena', 256)
Compute spherical coordinates $(\th,\phi)$.
v = X0 - repmat(mean(X0, 2), [1 n])
theta = acos(v(1, : )./ sqrt(sum(v.^2)))/ pi
phi = (atan2(v(2, : ), v(3, : ))/ pi + 1)/ 2
Interpolate the texture on the mesh to obtain $f_0$.
x = linspace(0, 1, size(M, 1))
f0 = interp2(x, x, M', theta, phi)'
Display the initial condition $f_0$ as a function on the mesh.
options.face_vertex_color = f0(: )
plot_mesh(X0, F, options)
lighting none
The PDE is discretized as $$ f^{(\ell+1)} = f^{(\ell)} - \tau \tilde L f^{(\ell)}. $$ Here $f^{(\ell)}$ is intended to approximate $f(\cdot,t)$ at time $t=\ell\tau$.
Final time.
Tmax = 50
The time step $\tau$ should be smaller than 1 to ensure stability of the scheme.
tau = .4
Number of iterations of the method.
niter = ceil(Tmax/ tau)
Initial solution at time t=0.
f = f0
One step of discretization.
f = f - tau*tL*f
Exercise 2
Compute the linear heat diffusion.
solutions.exo2()
## Insert your code here.
It is possible to apply the heat flow to each coordinate of the vertices $X$.
Initial solution at time t=0.
X = X0
We use an explicit discretization in time of the PDE. Here is one iteration. Note that since the matrix |X| stores the coordinates as row vectors, the application of $\tilde L$ is obtained by multiplication on the right by |tL'|.
X = X - tau*X*tL'
Exercise 3
Compute the linear heat diffusion by iterating this gradient descent.
solutions.exo3()
## Insert your code here.
The wave equation is the prototype of hyperbolic PDE. It produces a propagation of the initial data according to some speed. $$ \pdd{f}{t} = -L f \qandq \pd{f}{t}(\cdot,t=0)=g_0 \qandq f(\cdot,t=0)=f_0. $$
Compute initial conditions $f_0$: little Gaussian bumps as far as possible one from each other.
n = size(X0, 2)
width of the Gaussian
sigma = .002
numer of points
npoints = 16
initial condition
f0 = zeros(n, 1)
previous distances
D = zeros(n, 1) + 1e10
for i in 1: npoints:
% select a point farthest away
[tmp, p] = max(D)
% compute the distance beetween the point and the other vertices
d = sum((X0 - X0(: , p)*ones(1, n)).^2)
f0 = f0 + (-1)^mod(i, 2) *exp(-d(: )/ (2*sigma.^2))
D = min(D, d(: ))
Display it.
options.face_vertex_color = f0
plot_mesh(X0, F, options)
colormap(jet(256))
Set up parameters for the wave equation.
Tmax = 80
tau = .5
niter = round(Tmax/ tau)
We use explicite time stepping to evaluate the derivatives, which defines $$ f^{(\ell+1)} = 2 f^{(\ell)} - f^{(\ell-1)} - \tau^2 \tilde L f^{(\ell)}. $$
During the iteration, you need to keep track of the value |f1| of the solution of the PDE at the previous iteration. Since we assume that the wave startes with zero initial velocity $g_0=0$, we initialize as
f1 = f0
Perform a descent step of the equation.
update = lambda f, f1: deal(2*f - f1 - tau^2 * tL*f, f)
[f, f1] = update(f, f1)
Exercise 4
Solve the wave equation PDE.
solutions.exo4()
## Insert your code here.