# 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 `ChebyshevInterval`

s

`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.*