Poisson Equation

\[ \frac{\partial^2 u}{\partial x^2} = b \qquad x \in [0,1] \]

We solve only interior points: the endpoints are set to zero.

\[ u(0) = u(1) = 0, \qquad b = \sin(2\pi x) \]

using Plots, BenchmarkTools

function plot_solution(x, u)
    plot([0;x;1],[0;u;0], label="computed")
    scatter!([0;x;1],-sin.(2π*[0;x;1])/(4π^2),label="exact")
end

plot_solution (generic function with 1 method)

Δx = 0.05
x = Δx:Δx:1-Δx 
N = length(x)

A = zeros(N,N)
for i in 1:N, j in 1:N
    abs(i-j) <= 1 && (A[i,j] +=1)
    i==j          && (A[i,j] -=3)
end

B = sin.(2π*x) * Δx^2
u = A \ B

19-element Vector{Float64}:
 -0.007892189393343806
 -0.015011836300750243
 -0.020662020077425496
 -0.024289661368163386
 -0.025539661368163387
 -0.024289661368163386
 -0.0206620200774255
 -0.015011836300750247
 -0.007892189393343808
 -1.5770213417970973e-18
  0.007892189393343803
  0.015011836300750236
  0.02066202007742549
  0.024289661368163375
  0.025539661368163376
  0.02428966136816338
  0.020662020077425493
  0.015011836300750241
  0.007892189393343805

plot_solution(x, u)

SparseArrays

using SparseArrays

Δx = 0.05
x = Δx:Δx:1-Δx 
N = length(x)
B = sin.(2π*x) * Δx^2

19-element Vector{Float64}:
  0.0007725424859373687
  0.001469463130731183
  0.002022542485937369
  0.0023776412907378845
  0.0025000000000000005
  0.0023776412907378845
  0.002022542485937369
  0.0014694631307311833
  0.0007725424859373689
  3.0616169978683835e-19
 -0.0007725424859373683
 -0.0014694631307311829
 -0.002022542485937369
 -0.0023776412907378845
 -0.0025000000000000005
 -0.0023776412907378845
 -0.0020225424859373693
 -0.0014694631307311838
 -0.0007725424859373692

P = spdiagm( -1 =>    ones(Float64,N-1),
              0 => -2*ones(Float64,N),
              1 =>    ones(Float64,N-1))

19×19 SparseMatrixCSC{Float64, Int64} with 55 stored entries:
⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠈⠻⠆⎦

u = P \ B

19-element Vector{Float64}:
 -0.007892189393343805
 -0.015011836300750241
 -0.020662020077425496
 -0.024289661368163382
 -0.025539661368163383
 -0.024289661368163386
 -0.0206620200774255
 -0.015011836300750248
 -0.007892189393343811
 -6.308085367188389e-18
  0.0078921893933438
  0.01501183630075024
  0.020662020077425496
  0.024289661368163382
  0.025539661368163387
  0.024289661368163386
  0.020662020077425496
  0.015011836300750243
  0.007892189393343806

plot_solution(x, u)

LinearAlgebra

using LinearAlgebra

Δx = 0.05
x = Δx:Δx:1-Δx 
N = length(x)
B = sin.(2π*x) * Δx^2

19-element Vector{Float64}:
  0.0007725424859373687
  0.001469463130731183
  0.002022542485937369
  0.0023776412907378845
  0.0025000000000000005
  0.0023776412907378845
  0.002022542485937369
  0.0014694631307311833
  0.0007725424859373689
  3.0616169978683835e-19
 -0.0007725424859373683
 -0.0014694631307311829
 -0.002022542485937369
 -0.0023776412907378845
 -0.0025000000000000005
 -0.0023776412907378845
 -0.0020225424859373693
 -0.0014694631307311838
 -0.0007725424859373692

DU =  ones(Float64, N-1)
D = -2ones(Float64, N)
DL =  ones(Float64, N-1)

18-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

T = Tridiagonal(DL, D, DU)

19×19 Tridiagonal{Float64, Vector{Float64}}:
 -2.0   1.0    ⋅     ⋅     ⋅     ⋅   …    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
  1.0  -2.0   1.0    ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅    1.0  -2.0   1.0    ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅    1.0  -2.0   1.0    ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅    1.0  -2.0   1.0       ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅    1.0  -2.0  …    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅    1.0       ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   …    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅       1.0    ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      -2.0   1.0    ⋅     ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅       1.0  -2.0   1.0    ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   …    ⋅    1.0  -2.0   1.0    ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅    1.0  -2.0   1.0    ⋅ 
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅    1.0  -2.0   1.0
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅        ⋅     ⋅     ⋅     ⋅    1.0  -2.0

u = T \ B

19-element Vector{Float64}:
 -0.007892189393343806
 -0.015011836300750243
 -0.020662020077425496
 -0.024289661368163386
 -0.025539661368163387
 -0.024289661368163386
 -0.0206620200774255
 -0.015011836300750247
 -0.007892189393343808
 -1.5770213417970973e-18
  0.007892189393343803
  0.015011836300750236
  0.02066202007742549
  0.024289661368163375
  0.025539661368163376
  0.02428966136816338
  0.020662020077425493
  0.015011836300750241
  0.007892189393343805

plot_solution(x, u)

LAPACK

using LinearAlgebra

Δx = 0.05
x = Δx:Δx:1-Δx ## Solve only interior points: the endpoints are set to zero.
N = length(x)
B = sin.(2π*x) * Δx^2

19-element Vector{Float64}:
  0.0007725424859373687
  0.001469463130731183
  0.002022542485937369
  0.0023776412907378845
  0.0025000000000000005
  0.0023776412907378845
  0.002022542485937369
  0.0014694631307311833
  0.0007725424859373689
  3.0616169978683835e-19
 -0.0007725424859373683
 -0.0014694631307311829
 -0.002022542485937369
 -0.0023776412907378845
 -0.0025000000000000005
 -0.0023776412907378845
 -0.0020225424859373693
 -0.0014694631307311838
 -0.0007725424859373692

DU =   ones(Float64, N-1)
D  = -2ones(Float64, N)
DL =   ones(Float64, N-1)

18-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

LAPACK.gtsv!(DL, D, DU, B)

19-element Vector{Float64}:
 -0.007892189393343806
 -0.015011836300750243
 -0.020662020077425496
 -0.024289661368163386
 -0.025539661368163387
 -0.024289661368163386
 -0.0206620200774255
 -0.015011836300750247
 -0.007892189393343808
 -1.5770213417970973e-18
  0.007892189393343803
  0.015011836300750236
  0.02066202007742549
  0.024289661368163375
  0.025539661368163376
  0.02428966136816338
  0.020662020077425493
  0.015011836300750241
  0.007892189393343805

plot_solution(x, B)