Solving PDEs

Poisson equation

We solve the Poisson equation

\[Δu(x,y) = f(x,y)\]

on the rectangle -1..1 × -1..1, with zero boundary conditions:

using ApproxFun
using LinearAlgebra

d = (-1..1)^2
x,y = Fun(d)
f = exp.(-10(x+0.3)^2-20(y-0.2)^2)  # use broadcasting as exp(f) not implemented in 2D
A = [Dirichlet(d); Laplacian()]
u = A \ [zeros(∂(d)); f];

Using a QR Factorization reduces the cost of subsequent calls substantially

QR = qr(A)
u = QR \ [zeros(∂(d)); f];


import Plots
Plots.plot(u; st=:surface, legend=false)

Many PDEs have weak singularities at the corners, in which case it is beneficial to specify a tolerance to reduce the time, as:

\(A, [zeros(∂(d)); f]; tolerance=1E-6);

Helmholtz Equation

We solve

\[(Δ + 100)u(x,y) = 0\]

on the rectangle -1..1 × -1..1, subject to the condition that $u=1$ on the boundary of the domain.

using ApproxFun
using LinearAlgebra

The rectangular domain may be expressed as a tensor product of ChebyshevIntervals

d = ChebyshevInterval()^2;

We construct the operator that collates the differential operator and boundary conditions

L = [Dirichlet(d); Laplacian()+100I];

We compute the QR decomposition of the operator to speed up the solution

Q = qr(L);
ApproxFunBase.resizedata!(Q,:,4000);

The boundary condition is a function that is equal to one on each edge

boundary_cond = ones(∂(d));

Solve the differential equation. This function has weak singularities at the corner, so we specify a lower tolerance to avoid resolving these singularities completely

u = \(Q, [boundary_cond; 0.]; tolerance=1E-7);

Plot the solution

import Plots
Plots.plot(u; st=:surface, legend=false)

Heat equation

We solve the heat equation

\[\frac{\partial}{\partial t} u(x,t) = \frac{\partial^2}{\partial x^2} u(x,t)\]

on the spatial domain -1..1 from $t=0$ to $t=1$, with zero boundary conditions, and the initial condition $u(x,0) = u_0(x)$, where we choose $u_0(x)$ to be a sharp Gaussian.

using ApproxFun
using LinearAlgebra

dx = -1..1;
dt = 0..1;

We construct a 2D domain that is the tensor product of x and t

d = dx × dt;

The initial condition

u0 = Fun(x->exp(-x^2/(2*0.2^2)), dx);

We define the derivatives on the 2D domain

Dx = Derivative(d,[1,0]);
Dt = Derivative(d,[0,1]);

The operator along with the initial and the boundary conditions

A = [Dt - Dx^2; I⊗ldirichlet(dt); bvp(dx)⊗I];

We collect the initial condition along with zeros for the two boundary conditions and the equation

rhs = [0.0; u0; 0.0; 0.0];

Solve the equation with an arbitrary tolerance. A lower tolerance will increase accuracy at the expense of execution time

u = \(A, rhs, tolerance=1e-4);

Plot the solution

import Plots
xplot = -1:0.02:1
p = Plots.plot(xplot, u.(xplot, 0), label="t=0", legend=true, linewidth=2)
for t in [0.05, 0.1, 0.2, 0.5, 0.8]
	Plots.plot!(xplot, u.(xplot, t), label="t=$t")
end
p

This page was generated using Literate.jl.