xxxxxxxxxx
 
begin 
    using Pkg
    Pkg.activate(mktempdir())
    Pkg.add("PyPlot")
    Pkg.add("PlutoUI")
    Pkg.add("DataFrames")
    using PlutoUI
    using PyPlot
    using DataFrames
end

15.0 s

Matrices from partial differential equations

As we focus in this course on partial differential equations, we need discuss matrices which evolve from the discretization of PDEs.

Are there any structural or numerical patterns in these matrices we can take advantage of with regard to memory and time complexity when solving linear systems ?

In this lecture we introduce a relatively simple "drosophila" problem which we will use do discuss these issues.

For the start we use simple structured disceretization grids and a finite difference approach to the discretization. Later, this will be generalized to more general grids and to finite element and finite volume discretization methods.

2.0 ms

Heat conduction in a one-dimensional rod

Heat source $f (x)$
$v_{L}, v_{R}$ : ambient temperatures
$α$ : boundary heat transfer coefficient
Second order boundary value problem in $Ω = [0, 1]$ :

${\begin{cases} - u^{″} (x) & = f (x) in Ω \\ - u^{'} (0) + α (u (0) - v_{L}) & = 0 \\ u^{'} (1) + α (u (1) - v_{R}) & = 0 \end{cases}$

The solution $u$ describes the equilibrium temperature distribution. Behind the second derivative is Fouriers law and the continuity equation
In math, the boundary conditions are called "Robin" or "third kind". They describe a heat in/outflux proportional to the difference between rod end temperature and ambient temperature
Fix a number of discretization points N
Let $h = \frac{1}{N - 1}$
Let $x_{i} = (i - 1) h$ $i = 1 \dots N$ be discretization points

12.3 μs

xxxxxxxxxx
 
N=11

1.1 μs

plotgrid (generic function with 1 method)

107 μs

xxxxxxxxxx
 
plotgrid(N)

633 ms

Finite difference approximation

We can approximate continuous functions $f$ by piecewise linear functions defined by the values $f_{i} = f (x_{i})$ . Using more points yields a better approximation:

7.8 μs

xxxxxxxxxx
 
plotgrid(N,func=x->0.5*sin(8*x)^2)

196 ms

Let $u_{i}$ approximations for $u (x_{i})$ and $f_{i} = f (x_{i})$
We can use a finite difference approximation to approximate $u^{'} (x_{i + \frac{1}{2}}) \approx \frac{u_{i + 1} - u_{i}}{h}$
Same approach for second derivative: $u^{″} (x_{i}) = \frac{u^{'} (x_{i + \frac{1}{2}}) - u^{'} (x_{i - \frac{1}{2}})}{h}$
Finite difference approximation of the PDE:

$\begin{aligned} - u^{'} (0) + α (u (0) - v_{L}) & \approx \frac{1}{2 h} (u_{0} - u_{2}) + α (u_{1} - v_{L}) = 0 \\ - u^{″} (x_{i}) - f (x_{i}) & \approx \frac{- u_{i + 1} + 2 u_{i} - u_{i - 1}}{h^{2}} - f_{i} = 0 & (i = 1 \dots N) (*) \\ u^{'} (1) + α (u (1) - v_{R}) & \approx \frac{1}{2 h} (u_{N + 1} - u_{N - 1}) + α (u_{N} - v_{R}) = 0 \end{aligned}$

9.3 μs

xxxxxxxxxx
 
plotgrid(N,mirror=true)

114 ms

Here, we introduced "mirror values" $u_{0}$ and $u_{N + 1}$ in order to approximate the boundary conditions accurately, such that the finite difference formulas used to approximate $u^{'} (0)$ or $u^{'} (1)$ are centered around these values.
After rearranging, these values can be expressed via the boundary conditions:

$\begin{aligned} u_{0} & = u_{2} + 2 h α (u_{1} - v_{L}) \\ u_{N + 1} & = u_{N - 1} + 2 h α (u_{N} - v_{L}) \end{aligned}$

Finally, they can be replaced in $(*)$

5.9 μs

Then, the system after multiplying by $h$ is reduced to:

$\begin{aligned} \frac{1}{h} ((1 + h α) u_{1} - u_{2}) = & \frac{h}{2} f_{1} + α v_{L} \\ \frac{1}{h} (- u_{i + 1} + 2 u_{i} - u_{i - 1}) & = h f_{i} & (i = 2 \dots N - 1) \\ \frac{1}{h} ((1 + h α) u_{N} - u_{N - 1}) & = \frac{h}{2} f_{N} + α v_{R} \end{aligned}$

2.4 μs

The resulting discretization matrix is

$A = (\begin{matrix} α + \frac{1}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ ⋱ & ⋱ & ⋱ & ⋱ \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{1}{h} + α \end{matrix})$

Outside of the three diagonals, the entries are zero.

3.2 μs

The right hand side is:

$(\begin{matrix} \frac{h}{2} f_{1} + α v_{L} \\ h f_{2} \\ h f_{3} \\ ⋮ \\ h f_{N - 2} \\ h f_{N - 1} \\ \frac{h}{2} f_{N} + α v_{R} \end{matrix})$

2.5 μs

Let us define functions assembling these:

2.3 μs

heatmatrix1d (generic function with 1 method)

xxxxxxxxxx
 
function heatmatrix1d(N;α=1)
    A=zeros(N,N)
    h=1/(N-1)
    A[1,1]=1/h+α
    for i=2:N-1
        A[i,i]=2/h  
    end
    for i=1:N-1
        A[i,i+1]=-1/h
    end
    for i=2:N
        A[i,i-1]=-1/h
    end
    A[N,N]=1/h+α
    A
end

54.7 μs

heatrhs1d (generic function with 1 method)

xxxxxxxxxx
 
function heatrhs1d(N;vl=0,vr=0,func=x->0,α=1)
    h=1/(N-1)
    F=zeros(N)
    F[1]=h/2*func(0)+α*vl
    for i=2:N-1
        F[i]=h*func((i-1)*h)
    end
    F[N]=h/2*func(1)+α*vr
    F
end

123 μs

xxxxxxxxxx
 
α=100

959 ns

xxxxxxxxxx
 
N1=10000

1.2 μs

10000×10000 Array{Float64,2}:
 10099.0  -9999.0      0.0      0.0      0.0  …      0.0      0.0      0.0      0.0
 -9999.0  19998.0  -9999.0      0.0      0.0         0.0      0.0      0.0      0.0
     0.0  -9999.0  19998.0  -9999.0      0.0         0.0      0.0      0.0      0.0
     0.0      0.0  -9999.0  19998.0  -9999.0         0.0      0.0      0.0      0.0
     0.0      0.0      0.0  -9999.0  19998.0         0.0      0.0      0.0      0.0
     0.0      0.0      0.0      0.0  -9999.0  …      0.0      0.0      0.0      0.0
     0.0      0.0      0.0      0.0      0.0         0.0      0.0      0.0      0.0
     ⋮                                        ⋱                             
     0.0      0.0      0.0      0.0      0.0         0.0      0.0      0.0      0.0
     0.0      0.0      0.0      0.0      0.0  …  -9999.0      0.0      0.0      0.0
     0.0      0.0      0.0      0.0      0.0     19998.0  -9999.0      0.0      0.0
     0.0      0.0      0.0      0.0      0.0     -9999.0  19998.0  -9999.0      0.0
     0.0      0.0      0.0      0.0      0.0         0.0  -9999.0  19998.0  -9999.0
     0.0      0.0      0.0      0.0      0.0         0.0      0.0  -9999.0  10099.0

xxxxxxxxxx
 
A=heatmatrix1d(N1,α=α)

92.2 ms

Float641

5.0005e-5

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

0.00010001

more9991

0.00010001

9992

0.00010001

9993

0.00010001

9994

0.00010001

9995

0.00010001

9996

0.00010001

9997

0.00010001

9998

0.00010001

9999

0.00010001

10000

5.0005e-5

xxxxxxxxxx
 
b=heatrhs1d(N1,func=x->1,α=α)

25.4 ms

Float641

0.005

0.00505

0.00509999

0.00514997

0.00519994

0.0052499

0.00529985

0.00534979

0.00539972

0.00544964

0.00549955

0.00554945

0.00559934

0.00564922

0.00569909

0.00574895

0.0057988

0.00584864

0.00589847

0.00594829

more9991

0.00544964

9992

0.00539972

9993

0.00534979

9994

0.00529985

9995

0.0052499

9996

0.00519994

9997

0.00514997

9998

0.00509999

9999

0.00505

10000

0.005

xxxxxxxxxx
 
u=A\b

12.5 s

xxxxxxxxxx
 
begin
    clf()
    plot(collect(0:1/(N1-1):1),u)
    grid()
    gcf()
end

716 ms

For this example, we created an $N \times N$ matrix where all entries outside of the main diagonal and the two adjacent ones are zero:

Fraction of nonzero entries: 0.00029998
Ratio of nonzero entries to number of unknowns: 2.9998
In fact, this matrix has $O (N)$ nonzero entries.

331 ms

2D heat conduction

2.4 μs

Just pose the heat problem in a 2D domain $Ω = (0, 1) \times (0, 1)$ :

${\begin{cases} - \frac{\partial^{2} u}{\partial x^{2}} - \frac{\partial^{2} u}{\partial y^{2}} & = f (x, y) in Ω \\ \frac{\partial u}{\partial n} + α (u - v) & = 0 on \partial Ω \end{cases}$

We use 2D regular discretization $n \times n$ grid with grid points $x_{i j} = ((i - 1) h, (j - 1) h)$ . The finite difference approximation yields:

$\frac{- u_{i - 1, j} - u_{i, j - 1} + 4 u_{i j} - u_{i + 1, j} - u_{i, j + 1}}{h^{2}} = f_{i j}$

This just comes from summing up the 1D finite difference formula for the $x$ and $y$ directions.

We do not discuss the boundary conditions here.

The $n \times n$ grid leads to an $n^{2} \times n^{2}$ matrix!

4.6 μs

plotgrid2d (generic function with 1 method)

xxxxxxxxxx
 
function plotgrid2d(N;text=true, func=nothing)
    clf()
    ax=PyPlot.axes(aspect=1)
    x=[(i-1)/(N-1) for i=1:N]
    y=[(i-1)/(N-1) for i=1:N]
    
    for i=1:N
        plot([x[i],x[i]],[0,1],linewidth=1,color="k")
        plot([0,1],[y[i],y[i]],linewidth=1,color="k")
    end
    if func!=nothing
        f=[func(x[i],y[j]) for i=1:N, j=1:N]
        contourf(x,y,f,cmap="hot")
    end
    if text
      ij=1
      for j=1:N
        for i=1:N
          ax.text(x[i],y[j]-0.035,"\$x_{$(ij)}\$",fontsize=10,color=:blue)
            ij=ij+1
        end
      end
    end
    fig=PyPlot.gcf()
    fig.set_size_inches(5,5)
    fig
end

108 μs

xxxxxxxxxx
 
plotgrid2d(5)

160 ms

Matrix and right hand side assembly inspired by the finite volume method which will be covered later in the course. The result is the same as for the finite difference method with the mirror trick for the boundary condition.

2.3 μs

heatmatrix2d (generic function with 1 method)

xxxxxxxxxx
 
function heatmatrix2d(n;α=1)
    function update_pair(A,v,i,j)
        A[i,j]+=-v
        A[j,i]+=-v
        A[i,i]+=v
        A[j,j]+=v
    end
    N=n^2
    h=1.0/(n-1)
    A=zeros(N,N)
    l=1
        for j=1:n
            for i=1:n
                if i<n
                    update_pair(A,1.0,l,l+1)
                end
                if i==1|| i==n
                    A[l,l]+=α
                end
                if j<n
                    update_pair(A,1,l,l+n)
                end
                if j==1|| j==n
                    A[l,l]+=α
                end
                l=l+1
            end
        end
    A
end
​

65.3 μs

heatrhs2d (generic function with 1 method)

xxxxxxxxxx
 
function heatrhs2d(n; rhs=(x,y)->0,bc=(x,y)->0,α=1.0) 
    h=1.0/(n-1)
    x=collect(0:h:1)
    y=collect(0:h:1)
    N=n^2
    f=zeros(N)
    for i=1:n-1
        for j=1:n-1
            ij=(j-1)*n+i
            f[ij]+=h^2/4*rhs(x[i],y[j])
            f[ij+1]+=h^2/4*rhs(x[i+1],y[j])
            f[ij+n]+=h^2/4*rhs(x[i],y[j+1])
            f[ij+n+1]+=h^2/4*rhs(x[i+1],y[j+1])
        end
    end
    
    for i=1:n
        ij=i
        fac=h
        if i==1 || i==n 
            fac=h/2
        end
        f[ij]+=fac*α*bc(x[i],0)
        ij=i+(n-1)*n
        f[ij]+=fac*α*bc(x[i],1)
    end
    for j=1:n
        fac=h
        if j==1 || j==n 
            fac=h/2
        end
        ij=1+(j-1)*n
        f[ij]+=fac*α*bc(0,y[j])
        ij=n+(j-1)*n
        f[ij]+=fac*α*bc(1,y[j])
    end
    f
end

97.3 μs

xxxxxxxxxx
 
n=5

1.4 μs

Float641

0.0

0.0

0.0

0.0

0.0

0.0

0.03125

-0.0441942

0.03125

8.11834e-18

0.0

-0.0441942

0.0625

-0.0441942

-1.14811e-17

0.0

0.03125

-0.0441942

0.03125

8.11834e-18

0.0

8.11834e-18

-1.14811e-17

8.11834e-18

2.10904e-33

xxxxxxxxxx
 
b2=heatrhs2d(n,rhs=(x,y)->sin(3*π*x)*sin(3*π*y),α=α)

107 ms

25×25 Array{Float64,2}:
 202.0   -1.0    0.0    0.0    0.0   -1.0   0.0  …    0.0    0.0    0.0    0.0    0.0
  -1.0  103.0   -1.0    0.0    0.0    0.0  -1.0       0.0    0.0    0.0    0.0    0.0
   0.0   -1.0  103.0   -1.0    0.0    0.0   0.0       0.0    0.0    0.0    0.0    0.0
   0.0    0.0   -1.0  103.0   -1.0    0.0   0.0       0.0    0.0    0.0    0.0    0.0
   0.0    0.0    0.0   -1.0  202.0    0.0   0.0       0.0    0.0    0.0    0.0    0.0
  -1.0    0.0    0.0    0.0    0.0  103.0  -1.0  …    0.0    0.0    0.0    0.0    0.0
   0.0   -1.0    0.0    0.0    0.0   -1.0   4.0       0.0    0.0    0.0    0.0    0.0
   ⋮                                  ⋮          ⋱    ⋮                         
   0.0    0.0    0.0    0.0    0.0    0.0   0.0       0.0    0.0    0.0    0.0   -1.0
   0.0    0.0    0.0    0.0    0.0    0.0   0.0  …  202.0   -1.0    0.0    0.0    0.0
   0.0    0.0    0.0    0.0    0.0    0.0   0.0      -1.0  103.0   -1.0    0.0    0.0
   0.0    0.0    0.0    0.0    0.0    0.0   0.0       0.0   -1.0  103.0   -1.0    0.0
   0.0    0.0    0.0    0.0    0.0    0.0   0.0       0.0    0.0   -1.0  103.0   -1.0
   0.0    0.0    0.0    0.0    0.0    0.0   0.0       0.0    0.0    0.0   -1.0  202.0

xxxxxxxxxx
 
A2=heatmatrix2d(n,α=α)

38.0 ms

In order to inspect the matrix, we can turn it into a DataFrame, which can be browsed.

2.4 μs

	x1	x2	x3	x4	x5	x6	x7	x8
	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	202.0	-1.0	0.0	0.0	0.0	-1.0	0.0	0.0
2	-1.0	103.0	-1.0	0.0	0.0	0.0	-1.0	0.0
3	0.0	-1.0	103.0	-1.0	0.0	0.0	0.0	-1.0
4	0.0	0.0	-1.0	103.0	-1.0	0.0	0.0	0.0
5	0.0	0.0	0.0	-1.0	202.0	0.0	0.0	0.0
6	-1.0	0.0	0.0	0.0	0.0	103.0	-1.0	0.0
7	0.0	-1.0	0.0	0.0	0.0	-1.0	4.0	-1.0
8	0.0	0.0	-1.0	0.0	0.0	0.0	-1.0	4.0
9	0.0	0.0	0.0	-1.0	0.0	0.0	0.0	-1.0
10	0.0	0.0	0.0	0.0	-1.0	0.0	0.0	0.0
more
25	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

xxxxxxxxxx
 
DataFrame(A2)

317 ms

Float641

4.35673e-7

4.4003e-5

-6.20691e-5

4.4003e-5

4.35673e-7

4.4003e-5

0.00459394

-0.00648112

0.00459394

4.4003e-5

-6.20691e-5

-0.00648112

0.00914388

-0.00648112

-6.20691e-5

4.4003e-5

0.00459394

-0.00648112

0.00459394

4.4003e-5

4.35673e-7

4.4003e-5

-6.20691e-5

4.4003e-5

4.35673e-7

xxxxxxxxxx
 
u2=A2\b2

18.6 μs

xxxxxxxxxx
 
begin
    clf()
    h=1.0/(n-1)
    x=collect(0:h:1)
    y=collect(0:h:1)
    
    contourf(x,y,reshape(u2,n,n),cmap="hot")
    fig=gcf()
    fig.set_size_inches(5,5)
    fig
end

152 ms

In order to achieve this, we stored a matrix which has only five nonzero diagonals as a full $N \times N$ matrix, where $N = n^{2}$ :

Fraction of nonzero entries: 0.168
Ratio of nonzero entries to number of unknowns: 4.2
In fact, this matrix has $O (N)$ nonzero entries.

... there must be a better way!

32.1 μs