This Jupyter notebook serves as supplementary material to the Julia code from the book Numerical Methods for Scientific Computing. By and large, the snippets are verbatim from the book an occasional explicit output, plot statement, or variable declaration. These code snippets are minimal working toy algorithms meant to better understand the mathematics that goes into them. They are tools for tinkering and learning. Play with them and have fun. And, perhaps you can repurpose a few of them. The notebook is designed to be nonlinear—you can jump around. The LinearAlgebra.jl and Plots.jl packages are assumed to be imported.
using LinearAlgebra, Plots
bucket = "https://raw.githubusercontent.com/nmfsc/data/master/";
The notebook was written and tested using Julia 1.9.2. Several specialized packages are used throughout this notebook.
packages = ["LinearAlgebra","Plots","Images","DelimitedFiles","SparseArrays",
"Graphs","GraphPlot","FFTW","Primes","BenchmarkTools","DifferentialEquations",
"Optim","Dierckx","Interpolations","SpecialFunctions","Roots","OffsetArrays",
"Wavelets","Interact","LsqFit","GLM","DataFrames","Flux","FastGaussQuadrature",
"Sundials","SciMLOperators","CSV","ColorSchemes","JuMP","GLPK","MLDatasets",
"Arpack","ProgressMeter","IterativeSolvers","Polynomials","ParameterizedFunctions",
"ImageShow","NLsolve","DiffEqDevTools","LinRegOutliers","ModelingToolkit","NonlinearSolve"];
#import Pkg; Pkg.add(packages)
Part 1: Numerical Linear Algebra
Chapter 1: A Review of Linear Algebra
Chapter 2: Direct Methods for Linear Systems
Chapter 3: Inconsistent Systems
Chapter 4: Computing Eigenvalues
Chapter 6: Fast Fourier Transform
Part 2: Numerical Methods for Analysis
Chapter 7: Preliminaries
Chapter 8: Solutions to Nonlinear Equations
Chapter 9: Interpolation
Chapter 10: Approximating Functions
Chapter 11: Differentiation and Integration
Part 3: Numerical Differential Equations
Chapter 12: Ordinary Differential Equations
Chapter 13: Parabolic Equations
Chapter 16: Fourier Spectral Methods
Part 4: Solutions
Numerical Linear Algebra
Numerical Analysis
Numerical Differential Equations
hilbert(n) = [1/(i+j-1) for i=1:n, j=1:n]
using Images
[Gray.(1 .- abs.(hilbert(n)\hilbert(n))) for n ∈ (10,15,20,25,50)]
Gaussian elimination. The following function implements a naïve Gaussian elimination algorithm for a matrix A
and vector b
. We'll verify the code using a random matrix-vector pair.
function gaussian_elimination(A,b)
n = size(A,1)
for j in 1:n
A[j+1:n,j] /= A[j,j]
A[j+1:n,j+1:n] -= A[j+1:n,j:j].*A[j:j,j+1:n]
end
for i in 2:n
b[i:i] -= A[i:i,1:i-1]*b[1:i-1]
end
for i in n:-1:1
b[i:i] = ( b[i] .- A[i:i,i+1:n]*b[i+1:n] )/A[i,i]
end
return b
end
A = rand(8,8); b = rand(8,1)
[A\b gaussian_elimination(A,b)]
ϵ = 1e-15; A = [1 1 2;2 2+ϵ 0; 4 14 4]; b = [-5;10;0.0]
A\b, gaussian_elimination(A,b)
Simplex method. The following three functions (get_pivot
, row_reduce
, and simplex
) implement a naïve simplex method. Let's use them to solve the LP problem "Find the maximum of the objective function $2x + y + z$ subject to the constraints $2x+ z \leq 3$, $4x+y + 2z \leq 2$, and $x+y \leq 1$" along with the dual LP problem. To better illustrate the simplex method, it may be helpful to add the command display(round.(tableau,digits=2))
to the while
loop.
function row_reduce!(tableau)
(i,j) = get_pivot(tableau)
G = tableau[i:i,:]/tableau[i,j]
tableau[:,:] -= tableau[:,j:j]*G
tableau[i,:] = G
end
function get_pivot(tableau)
j = argmax(tableau[end,1:end-1])
a, b = tableau[1:end-1,j], tableau[1:end-1,end]
k = findall(a.>0)
i = k[argmin(b[k]./a[k])]
return(i,j)
end
function simplex(c,A,b)
(m,n) = size(A)
tableau = [[A I b] ; [c' zeros(1,m) 0]]
while (any(tableau[end,1:n].>0))
row_reduce!(tableau)
end
p = findall(tableau[end,1:n].==0)
x = zeros(n,1)
[x[i]=tableau[:,i]⋅tableau[:,end] for i∈p]
z = -tableau[end,end]
y = -tableau[end,n.+(1:m)]
return((z=z, x=x, y=y))
end
A = [2 0 1;4 1 2;1 1 0]; c = [2;1;1]; b = [3;2;1]
solution = simplex(c,A,b)
using DelimitedFiles, SparseArrays
function get_adjacency_matrix(filename)
ij = readdlm(download(bucket*filename*".csv"),',',Int)
sparse(ij[:,1],ij[:,2],one.(ij[:,1]))
end
using Graphs, GraphPlot
g = get_adjacency_matrix("dolphins") |> Graph
gplot(g,layout=spectral_layout)
gplot(g,layout=spring_layout)
gplot(g,layout=circular_layout)
Revised simplex method.
using SparseArrays
function revised_simplex(c,A,b)
(m,n) = size(A)
ξ = i -> (z=spzeros(m);z[i]=1;z)
N = Vector(1:n); B = Vector(n .+ (1:m))
A = [A sparse(I, m, m)]
ABinv = sparse(I, m, m)
c = [c;spzeros(m,1)]
while(true)
j = findfirst(x->x>0,(c[N]'.-(c[B]'*ABinv)*A[:,N])[:])
if isnothing(j); break; end
q = ABinv*A[:,N[j]]
k = findall(q.>0)
i = k[argmin(ABinv[k,:]*b./q[k])]
B[i], N[j] = N[j], B[i]
ABinv -= ((q - ξ(i))/q[i])*ABinv[i:i,:]
end
i = findall(B.≤n)
x = zeros(n,1)
x[B[i]] = ABinv[i,:]*b
return((z=c[1:n]'*x, x=x, y=c[B]'*ABinv))
end
A = [2 0 1;4 1 2;1 1 0]; c = [2;1;1]; b = [3;2;1]
solution = simplex(c,A,b)
Zipf's law. Let's use ordinary least squares to find Zipf's law coefficients for the canon of Sherlock Holmes.
using DelimitedFiles
T = readdlm(download(bucket*"sherlock.csv"), '\t')[:,2]
n = length(T)
A = [ones(n,1) log.(1:n)]
B = log.(T)
c = A\B
print("ordinary least squares:\n"*string(c)*"\n")
zipf(x,c) = exp(c[1])*x.^c[2]
scatter(T,xaxis=:log,yaxis=:log,ma=0.4,msw=0)
plot!([1,n],zipf([1,n],c),w=3)
Constrained least squares. The constrained least squares problem of solving $\mathbf{Ax} = \mathbf{b}$ with the constraint condition $\mathbf{Cx}=\mathbf{d}$:
function constrained_lstsq(A,b,C,d)
x = [A'*A C'; C zeros(size(C,1),size(C,1))]\[A'*b;d]
x[1:size(A,2)]
end
function tls(A,B)
n = size(A,2)
V = svd([A B]).V
-V[1:n,n+1:end]/V[n+1:end,n+1:end]
end
A = [2 4; 1 -1; 3 1; 4 -8]; b = [1; 1; 4; 1];
xₒₗₛ = A\b; xₜₗₛ = tls(A,b)
print("ordinary least squares:"*string(xₒₗₛ))
print("\ntotal least squares:"*string(xₜₗₛ))
Image compression. Let's use singular value decomposition to compress an image. We'll choose a nominal rank k = 20
for demonstration. We'll use the Frobenius norm to compute the total pixelwise error in the compressed image. Then, we'll plot out all the singular values for comparison.
using Images
A = load(download(bucket*"laura.png")) .|> Gray
U, σ, V = svd(A);
k = 20
Aₖ = U[:,1:k] * Diagonal(σ[1:k]) * V[:,1:k]' .|> Gray
[A Aₖ]
norm(A-Aₖ) ≈ norm(σ[k+1:end])
ϵ² = 1 .- cumsum(σ)/sum(σ); scatter(ϵ²,xaxis=:log)
Nonnegative matrix factorization. The following code is a naive implementation of nonnegative matrix factorization using multiplicative updates (without a stopping criterion):
function nmf(X,p=6)
W = rand(Float64, (size(X,1), p))
H = rand(Float64, (p,size(X,2)))
for i in 1:50
W = W.*(X*H')./(W*(H*H') .+ (W.≈0))
H = H.*(W'*X)./((W'*W)*H .+ (H.≈0))
end
(W,H)
end
using Images
A = Gray.(load(download(bucket*"laura.png")))
W,H = nmf(Float64.(A),20)
[A Gray.(W*H)]
Eigenvalue condition number. The function condeig
computes the eigenvalue condition number. Let's use it on a small random matrix.
function condeig(A)
Sᵣ = eigvecs(A)
Sₗ = inv(Sᵣ')
Sₗ ./= sqrt.(sum(abs.(Sₗ.^2), dims=1))
1 ./ abs.(sum(Sᵣ.*Sₗ, dims=1))
end
condeig(rand(4,4))
PageRank. The following minimal code computes the PageRank of the very small graph by using the power method over 9 iterations
H = [0 0 0 0 1; 1 0 0 0 0; 1 0 0 0 1; 1 0 1 0 0; 0 0 1 1 0]
v = all(H.==0,dims=1)
H = H ./ (sum(H,dims=1)+v)
n = size(H,1)
d = 0.85
x = ones(n,1)/n
for i in 1:9
x = d*(H*x) .+ d/n*(v*x) .+ (1-d)/n
end
x
Shuffling the FFT matrix. Shuffling the rows of the F₁₂₈ yields four matrices that identical to (or diagonally similar) to F₆₄. Shuffling the rows of the F₆₄ yields four matrices that identical to (or diagonally similar) to F₃₂. And so on.
using FFTW, Images, ColorSchemes
perfectshuffle= A -> [A[1:2:end,:];A[2:2:end,:]]
F = real.(fft(I(128),1)); PF = F
for i = 1:6
F = perfectshuffle(F); PF = [PF F]
end
get(colorschemes[:gray1], PF, :extrema)
Radix-2 FFT. This chapter introduces several naive functions. The radix-2 FFT algorithm is written as a recursive function fftx2
and the inverse FFT is written as ifftx2
.
function fftx2(c)
n = length(c)
ω = exp(-2im*π/n)
if mod(n,2) == 0
k = collect(0:n/2-1)
u = fftx2(c[1:2:n-1])
v = (ω.^k).*fftx2(c[2:2:n])
return([u+v; u-v])
else
k = collect(0:n-1)
F = ω.^(k*k')
return(F*c)
end
end
ifftx2(y) = conj(fftx2(conj(y)))/length(y);
Fast Toeplitz multiplication. The function fasttoeplitz
computes the Toeplitz multiplication by padding out a Toeplitz matrix with zeros to make it circulant.
function fasttoeplitz(c,r,x)
n = length(x)
Δ = nextpow(2,n) - n
x₁ = [c; zeros(Δ); r[end:-1:2]]
x₂ = [x; zeros(Δ+n-1)]
ifftx2(fftx2(x₁).*fftx2(x₂))[1:n]
end
using FFTW, Primes, Plots
N = 10000
smooth(n,N) = (1:N)[all.(x->x<=n,factor.(Set,1:N))]
t₁ = [(x = randn(n); (n,@elapsed(fft(x)))) for n∈primes(N)]
t₂ = [(x = randn(n); (n,@elapsed(fft(x)))) for n∈smooth(5,N)]
plot(t₁,label="prime"); plot!(t₂,label="5-smooth")
plot!(ylims=(0,t₁[end][2]*3))
function bluestein_fft(x)
n = length(x)
ω = exp.((1im*π/n)*(0:n-1).^2)
ω.*fasttoeplitz(conj(ω),conj(ω),ω.*x)
end
using Primes
function primitiveroot(n)
ϕ = n - 1
p = factor(Set,ϕ)
for r = 2:n-1
all( [powermod(r, ϕ÷pᵢ, n) for pᵢ∈p] .!= 1 ) && return(r)
end
end
function rader_fft(x)
n = length(x)
r = primitiveroot(n)
P₊ = powermod.(r, 0:n-2, n)
P₋ = circshift(reverse(P₊),1)
ω = exp.((2im*π/n)*P₋)
c = x[1] .+ ifft(fft(ω).*fft(x[2:end][P₊]))
[sum(x); c[reverse(invperm(P₋))]]
end
Fast Poisson solver. The following code solves the Poisson equation using a naive fast Poisson solver and then compares the solution with the exact solution.
using FFTW
dst(x) = FFTW.r2r(x,FFTW.RODFT00)
idst(x) = dst(x)/(2^ndims(x)*prod(size(x).+1))
n = 50; x = (1:n)/(n+1); Δx = 1/(n+1)
v = 2 .- 2cos.(x*π)
λ = [v[i]+v[j]+v[k] for i∈1:n, j∈1:n, k∈1:n]./Δx^2
f = [(x-x^2)*(y-y^2) + (x-x^2)*(z-z^2)+(y-y^2)*(z-z^2)
for x∈x,y∈x,z∈x]
u = idst(dst(f)./λ);
norm(u - [(x-x^2)*(y-y^2)*(z-z^2)/2 for x∈x,y∈x,z∈x])
Image compression. Let's consider the following code that compresses a grayscale image A
to a factor of d
from its original size.
DCT image compression.
using Images, FFTW, SparseArrays
function dctcompress(A,d)
B = dct(Float64.(A))
idx = d*prod(size(A)) |> floor |> Int
tol = sort!(abs.(B[:]),rev=true)[idx]
B₀ = droptol!(sparse(B),tol)
idct(B₀) |> clamp01! .|> Gray
end
A = Gray.(load(download(bucket*"laura_square.png")))
[A dctcompress(A,0.05)]
bitstring(Float64(π))
Fast inverse square root. The following function implements the circa 1999 Q_rsqrt algorithm to approximate the reciprocal square root of a number.
function Q_rsqrt(x::Float32)
i = reinterpret(Int32,Float32(x))
i = Int32(0x5f3759df) - (i>>1)
y = reinterpret(Float32,i)
y = y * (1.5f0 - (0.5f0 * x * y * y))
end
Q_rsqrt(Float32(0.01))
using BenchmarkTools
@btime Q_rsqrt(Float32(x)) setup=(x=rand()) seconds=3
@btime 1/sqrt(x) setup=(x=rand()) seconds=3
Rump's catastrophic cancellation. The answer should be -0.827396...
What does Julia come up with?
a = 77617; b = 33096
333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b)
NaN
can be used to lift "pen off paper" when plotting a series of connected points.
We start with simple implementation of the bisection method.
function bisection(f,a,b,tolerance)
while abs(b-a) > tolerance
c = (a+b)/2
sign(f(c)) == sign(f(a)) ? a = c : b = c
end
(a+b)/2
end
bisection(x->sin(x),2,4,1e-14)
The Mandelbrot set. The following function takes the array bb
for the lower-left and upper-right corners of the bounding box, xn
for the number of horizontal pixels, and n
for the maximum number of iterations. The function returns a two-dimensional array M
that counts the number of iterations k
to escape $\{z\in\mathbb{C} \mid |z^{(k)}|>2\}$.
function escape(n, c=0, z=0)
for k in 1:n
z = z^2 + c
abs2(z) > 4 && return k
end
return n
end
function mandelbrot(bb, xn=800, n=200, z=0)
yn = (xn*(bb[4]-bb[2]))÷(bb[3]-bb[1]) |> Int
c = LinRange(bb[1],bb[3],xn)' .+ im*LinRange(bb[4],bb[2],yn)
[escape(n,c,z) for c∈c]
end
using Images
M = mandelbrot([-0.172, 1.0228, -0.1494, 1.0443])
save("mandelbrot.png",1 .- M./maximum(M))
load("mandelbrot.png")
Imaging the Julia set uses almost identical code. The Mandelbrot set lives in the $c$-domain with an initial value $z=0$, and the Julia set lives in the $z$-domain with a given value $c$. So the code for the Julia set requires only swapping the variables z
and c
.
function julia(bb, xn=800, n=200, c=-0.123+0.745im)
yn = (xn*(bb[4]-bb[2]))÷(bb[3]-bb[1]) |> Int
z = LinRange(bb[1],bb[3],xn)' .+ im*LinRange(bb[4],bb[2],yn)
[escape(n,c,z) for z∈z]
end
J = julia([-2, -2, 2, 2])
save("julia.png",1 .- J./maximum(J))
load("julia.png")
Homotopy continuation. The following snippet of code finds a root of $$x^3-3xy^2-1 =0$$ $$y^3-3x^2y = 0$$ with an initial guess $(x,y) = (1,1)$ using homotopy continuation:
using DifferentialEquations
f = z -> ((x,y)=tuple(z...); [x^3-3x*y^2-1; y^3-3x^2*y] )
df = (z,p,_) -> ((x,y)=tuple(z...);
-[3x^2-3y^2 -6x*y; -6x*y 3y^2-3x^2]\p)
z₀ = [1,1]
sol = solve(ODEProblem(df,z₀,(0,1),f(z₀)))
sol.u[end]
function plot_trace(x,f)
contour(-2:0.02:1, -0.25:0.02:3.25,(x,y) -> cbrt(f([x,y])), color=:red, colorbar=:none)
plot!(x[1,:],x[2,:],color=:black, linewidth=2, marker=:o, legend=:none)
end
∇f = x-> [-2(1-x[1])-400x[1]*(x[2]-x[1]^2); 200(x[2]-x[1]^2)]
x = [-1.8,3.0]; p = [0,0]; α = 0.001; β = 0.9
s = x
for i = 1:100
p = -∇f(x) + β*p
x += α*p
append!(s,x)
end
f = x -> (1-x[1])^2 + 100(x[2] - x[1]^2)^2
plot_trace(reshape(s, 2, :),f)
using Optim
x₀ = [-1.8,3.0]
options = Optim.Options(store_trace = true,extended_trace = true)
result = optimize(f, x₀, BFGS(),options)
s = Optim.x_trace(result); s = vcat(s...)
plot_trace(reshape(s, 2, :),f)
∇f = x -> [-2(1-x[1])-400x[1]*(x[2]-x[1]^2); 200(x[2]-x[1]^2)]
x = [-1.8,3.0]; α = 0.001; p = [0,0]; β = 0.9
for i = 1:100
p = -∇f(x) + β*p
x += α*p
end
using Optim
f = x -> (1-x[1])^2 + 100(x[2] - x[1]^2)^2
x₀ = [-1.8,3.0]
result = optimize(f, x₀, GradientDescent())
x = Optim.minimizer(result)
Splines. The function spline_natural
computes the coefficients m
of a cubic spline with natural boundary conditions through the nodes given by the arrays x
and y
. The function evaluate_spline
returns a set of n
points along the spline. The final code snippet tests these functions with several randomly selected points.
function spline_natural(xᵢ,yᵢ)
h = diff(xᵢ)
γ = 6*diff(diff(yᵢ)./h)
α = h[2:end-1]
β = 2(h[1:end-1]+h[2:end])
[0;SymTridiagonal(β,α)\γ;0]
end
function evaluate_spline(xᵢ,yᵢ,m,n)
h = diff(xᵢ)
B = yᵢ[1:end-1] .- m[1:end-1].*h.^2/6
A = diff(yᵢ)./h - h./6 .*diff(m)
x = range(minimum(xᵢ),maximum(xᵢ),length=n+1)
i = sum(xᵢ'.≤x,dims=2)
i[end] = length(xᵢ)-1
y = @. (m[i]*(xᵢ[i+1]-x)^3 + m[i+1]*(x-xᵢ[i])^3)/(6h[i]) +
A[i]*(x-xᵢ[i]) + B[i]
return(x,y)
end
x = LinRange(0,1,8); y = rand(8)
m = spline_natural(x,y)
X,Y = evaluate_spline(x,y,m,100)
scatter(x,y); plot!(X,Y)
The Julia Dierckx.jl package provides a wrapper of the dierckx Fortran library. This package provides an easy interface for building splines. Let's create a function and select some knots.
g = x-> @. max(1-abs(3-x),0)
xᵢ = 0:5; yᵢ = g(xᵢ)
x = LinRange(0,5,101);
using Dierckx, Plots
spline = Spline1D(xᵢ,yᵢ)
plot(x,spline(x)); plot!(x,g(x)); scatter!(xᵢ,yᵢ)
The Interpolations.jl package, which is still under development, can be used to evaluate up through third-order B-splines, but the nodes must be equally spaced. Let's fit a spline through the knots that we generated earlier.
using Interpolations
method = BSpline(Cubic(Natural(OnGrid())))
spline = scale(interpolate(yᵢ, method), xᵢ)
plot(x,spline(x)); plot!(x,g(x)); scatter!(xᵢ,yᵢ)
Bézier curves. The following function builds a Bernstein matrix. We'll then test the function on a set of points to create a cubic Bézier curve with a set of four randomly selected control points.
bernstein(n,t) = @. binomial(n,0:n)'*t^(0:n)'*(1-t)^(n:-1:0)'
n = 3
t = LinRange(0,1,100)
p = rand(n+1,2)
z = bernstein(n,t)*p
plot(p[:,1],p[:,2],marker=:o,opacity=0.3);
plot!(z[:,1],z[:,2],legend=:none)
n = 20; f = t -> sqrt(t)
y = bernstein(n,t)*f.(LinRange(0,1,n+1))
plot(t,y); plot!(t,f.(t))
Legendre polynomials. We can evaluate a Legendre polynomial of order $n$ using Bonnet's recursion formula.
function legendre(x,n)
n==0 && return(one.(x))
n==1 && return(x)
x.*legendre(x,n-1) .- (n-1)^2/(4(n-1)^2-1)*legendre(x,n-2)
end
x = LinRange(-1,1,100)
plot()
for n in 0:4
plot!(x,legendre(x,n))
end
current()
Chebyshev polynomials. We'll construct a Chebyshev differentiation matrix and use the matrix to solve a few simple problems.
function chebdiff(n)
x = -cos.(LinRange(0,π,n))
c = [2;ones(n-2);2].*(-1).^(0:n-1)
D = c./c'./(x .- x' + I)
D - Diagonal([sum(D,dims=2)...]), x
end
n = 15
D,x = chebdiff(n)
u = exp.(-4x.^2)
plot(x,D*u,marker=:o)
t = LinRange(-1,1,200)
plot!(t,-8t.*exp.(-4t.^2))
B = zeros(1,n); B[1] = 1
plot(x,[B;D]\[2;u],marker=:o)
Let's use the Chebyshev differentiation matrix to solve a boundary value problem (the Airy equation): $y''-256xy=0$ with $y(-1)=2$ and $y(1)=1$.
n = 15; k² = 256
D,x = chebdiff(n)
L = (D^2 - k²*Diagonal(x))
L[[1,n],:] .= 0; L[1,1] = L[n,n] = 1
y = L\[2;zeros(n-2);1]
plot(x,y,marker=:o)
The analytical solution can be expressed as the sum of Airy functions.
using SpecialFunctions
a = [airyai(-∛k²) airybi(-∛k²);airyai(∛k²) airybi(∛k²)]\[2;1]
airy(x) = a[1]*airyai(∛k²*x) + a[2]*airybi(∛k²*x)
scatter(x,y)
plot!(t,airy.(t))
Even with just 15 Chebyshev points, the numerical solution is a pretty good match to the analytical solution. Let's plot the \ell^\infty$-error as a function of the number of Chebyshev points. How many points do we need until we reach machine precision?
λ² = 16^2; N = 6:2:60; ϵ = []
for n ∈ N
D,x = chebdiff(n)
L = (D^2 - k²*Diagonal(x))
L[[1,n],:] .= 0; L[1,1] = L[n,n] = 1
y = L\[2;zeros(n-2);1]
append!(ϵ,norm(y - airy.(x),Inf))
end
plot(N,ϵ,yaxis=:log,marker=:o,legend=:none)
Wavelets. The function scaling
returns the scaling function (father wavelet). We can use it to generate the wavelet function (mother wavelet). As an example, we will plot the Daubechies $D_4$ with $c_k = (1+\sqrt{3},3+\sqrt{3},3-\sqrt{3},1-\sqrt{3}])/4$ and $\phi(k) = (0,1+\sqrt{3},1-\sqrt{3},0)/2$.
using OffsetArrays
function scaling(c,ϕₖ,n)
m = length(c)-1; ℓ = 2^n
ϕ = OffsetVector(zeros(3*m*ℓ),-m*ℓ)
k = (0:m)*ℓ
ϕ[k] = ϕₖ
for j = 1:n
for i = 1:m*2^(j-1)
x = (2i-1)*2^(n-j)
ϕ[x] = c ⋅ ϕ[2x .- k]
end
end
((0:m*ℓ-1)/ℓ, ϕ[0:m*ℓ-1])
end
c = [1+√3, 3+√3, 3-√3, 1-√3]/4
z = [0, 1+√3, 1-√3, 0]/2
(x,ϕ) = scaling(c, z, 8)
plot(x,ϕ)
ψ = zero(ϕ); n = length(c)-1; ℓ = length(ψ)÷2n
for k∈0:n
ψ[(k*ℓ+1):(k+n)*ℓ] += (-1)^k*c[n-k+1]*ϕ[1:2:end]
end
plot(x,ψ)
DWT image compression. The following code explores using a discrete wavelet transform along with filtering as means of image compression.
using Wavelets, Images
img = float.(Gray.(load(download(bucket*"laura_square.png"))))
B = dwt(Float32.(img), wavelet(WT.haar))
[img Gray.(1 .- clamp01.(sqrt.(abs.(B))))]
function dwtfilter(channels,wt,k)
map(1:3) do i
A = dwt(channels[i,:,:], wavelet(wt))
threshold!(A, HardTH(), k)
clamp01!(idwt(A, wavelet(wt)))
end
end
using Interact, Wavelets, Images
img = float.(load(download(bucket*"laura_square.png")))
channels = channelview(img)
func = Dict("Haar"=>WT.haar, "D₄"=>WT.db4, "Coiflet"=>WT.coif4)
@manipulate for wt in togglebuttons(func; label="Transform"),
k in slider(0:0.05:1.0; value = 0, label = "Threshold")
colorview(RGB,dwtfilter(channels,wt,k)...)
end
Nonlinear least squares approximation. The function gauss_newton
solves a nonlinear least squares problem using the Levenberg–Marquardt method. The Jacobian is approximated numerically by the function jacobian
using the complex-step method. The solver is then used to find the parameters for a two-Gaussian problem.
function jacobian(f,x,c)
J = zeros(length(x),length(c))
for k in (n = 1:length(c))
J[:,k] .= imag(f(x,c+1e-8im*(k .== n)))/1e-8
end
return J
end
function gauss_newton(f,x,y,c)
r = y - f(x,c)
for j = 1:100
G = jacobian(f,x,c)
M = G'*G
c += (M+Diagonal(M))\(G'*r)
norm(r-(r=y-f(x,c))) < 1e-12 && return(c)
end
print("Gauss-Newton did not converge.")
end
f = (x,c) -> @. c[1]*exp(-c[2]*(x-c[3])^2) +
c[4]*exp(-c[5]*(x-c[6])^2)
x = 8*rand(100)
y = f(x,[1, 3, 3, 2, 3, 6]) + 0.1*randn(100);
c₀ = [2, 0.3, 2, 1, 0.3, 7]
c = gauss_newton(f,x,y,c₀)
if !isnothing(c)
X = LinRange(0,8,200)
scatter(x,y,opacity=0.5)
plot!(X,f(X,c))
end
In practice, we might use the LsqFit.jl package which solves nonlinear least squares problems.
using LsqFit
cf = curve_fit(f,x,y,c₀)
c, r = cf.param, cf.resid
scatter(x,y,opacity=0.5)
plot!(X,f(X,c))
Logistic regression. We'll first define the logistic function and generate synthetic data. Then, we apply Newton's method. Finally, we plot the fit logistic function.
σ = x -> @. 1/(1+exp(-x))
x = rand(30); y = ( rand(30) .< σ(10x.-5) );
X, w = [one.(x) x], zeros(2,1)
for i=1:10
S = σ(X*w).*(1 .- σ(X*w))
w += (X'*(S.*X))\(X'*(y - σ(X*w)))
end
scatter(x,y)
t = LinRange(minimum(x),maximum(x),50)
plot!(t,σ([one.(t) t]*w))
In practice, we might use Julia's GLM library to solve the logistic regression problem.
using GLM, DataFrames
data = DataFrame(X=x, Y=y)
logit = glm(@formula(Y ~ X), data, Bernoulli(), LogitLink())
Neural Networks. Let's use a neural network to find a function that approximates semi-circular data.
N = 100; θ = LinRange(0,π,N)'
x = cos.(θ); x̃ = [one.(x);x]
y = sin.(θ) + 0.05*randn(1,N)
n = 20; W₁ = rand(n,2); W₂ = randn(1,n)
ϕ = x -> max.(x,0)
dϕ = x -> (x.>0)
α = 0.1
for epoch = 1:1000
ŷ = W₂ * ϕ(W₁*x̃)
∂L∂y = 2(ŷ-y)/N
∂L∂W₁ = dϕ(W₁*x̃) .* (W₂' * ∂L∂y) * x̃'
∂L∂W₂ = ∂L∂y * ϕ(W₁*x̃)'
W₁ -= 0.1α * ∂L∂W₁
W₂ -= α * ∂L∂W₂
end
scatter(x̃[2,:],y',opacity=0.3)
x̂ = LinRange(-1.2,1.2,200)'; x̂ = [one.(x̂);x̂]; ŷ = W₂ * ϕ(W₁*x̂)
plot!(x̂[2,:],ŷ',width=3)
We can solve the same problem using a sigmoid activation function:
n = 20; W₁ = rand(n,2); W₂ = randn(1,n)
ϕ = x -> @. 1 / (1 + exp(-x))
dϕ = x -> @. ϕ(x)*(1-ϕ(x))
α = 0.1
for epoch = 1:10000
ŷ = W₂ * ϕ(W₁*x̃)
L = norm(ŷ - y)
∂L∂y = (ŷ-y)/L
∂L∂W₁ = dϕ(W₁*x̃) .* (W₂' * ∂L∂y) * x̃'
∂L∂W₂ = ∂L∂y * ϕ(W₁*x̃)'
W₁ -= 0.1α * ∂L∂W₁
W₂ -= α * ∂L∂W₂
end
scatter(x̃[2,:],y',opacity=0.3)
x̂ = LinRange(-1.2,1.2,200)'; x̂ = [one.(x̂);x̂]; ŷ = W₂ * ϕ(W₁*x̂)
plot!(x̂[2,:],ŷ',width=3)
In practice, we might use Flux.jl to solve the same problem. The recipe is as before: build a model with initially unknown parameters, prescribe a loss function, add some data, choose an optimizer such as gradient descent, and then iterate.
using Flux
model = Chain(Dense(1, n, relu), Dense(n, 1))
loss(x,y) = Flux.Losses.mse(model(x), y)
parameters = Flux.params(model)
data = [(x,y)]
optimizer = Descent(0.1)
for epochs = 1:10000
Flux.train!(loss, parameters, data, optimizer)
end
scatter(x',y',ma=0.4,msw=0)
plot!(x',model(x)',width=3)
Let's compute the coefficients to the third-order approximation to $f'(x)$ using nodes at $x-h$, $x$, $x+h$ and $x+2h$.
δ = [-1,0,1,2]
n = length(δ)
V = δ.^(0:n-1)' .// factorial.(0:n-1)'
C = inv(V)
function richardson(f,x,m,n=m)
n > 0 ?
(4^n*richardson(f,x,m,n-1) - richardson(f,x,m-1,n-1))/(4^n-1) :
ϕ(f,x,2^m)
end
ϕ = (f,x,n) -> (f(x+1/n) - f(x-1/n))/(2/n)
richardson(x->sin(x),0,4)
struct Dual
value
deriv
end
Dual(x) = Dual(x,1)
value(x) = x isa Dual ? x.value : x
deriv(x) = x isa Dual ? x.deriv : 0
Base.:+(u, v) = Dual(value(u)+value(v), deriv(u)+deriv(v))
Base.:-(u, v) = Dual(value(u)-value(v), deriv(u)-deriv(v))
Base.:*(u, v) = Dual(value(u)*value(v),
value(u)*deriv(v)+value(v)*deriv(u))
Base.:sin(u) = Dual(sin(value(u)), cos(value(u))*deriv(u))
x = Dual(0)
y = x + sin(x)
Now, let's apply the code above to compute the Jacobian of the system $$y_1 = x_1x_2 + \sin x_2$$$$y_2 = x_1x_2 - \sin x_2$$ evaluated at $(x_1,x_2) = (2,\pi)$.
x₁ = Dual(2,[1 0])
x₂ = Dual(π,[0 1])
y₁ = x₁*x₂ + sin(x₂)
y₂ = x₁*x₂ - sin(x₂)
value(y₁),value(y₁),deriv(y₁),deriv(y₂)
Romberg method. We can use the following trapezoidal quadrature to make a Romberg method using Richardson extrapolation. We first define the function trapezoidal
for composite trapezoidal quadrature. By redefining phi
to be trapezoidal
, we can simply apply the function D
that we used to define Richardson extrapolation. We'll verify the code by integrating $\sin x$ from $0$ to $\pi/2$.
function trapezoidal(f,x,n)
F = f.(LinRange(x[1],x[2],n+1))
(F[1]/2 + sum(F[2:end-1]) + F[end]/2)*(x[2]-x[1])/n
end
ϕ = (f,x,n) -> trapezoidal(f,x,n)
richardson(x->sin(x),[0,π/2],4)
Composite trapezoidal method. Let's examine the convergence rate for the composite trapezoidal rule applied to the function $f(x) = x + (x - x^2)^p$ over the interval $[0,2]$ with $p = 1,2,\dots,7$. We can do this by finding the logl-og slope of the error as a function of subintervals $n$. We find that the error of the trapezoidal rule is $O(n^2)$ when $p=1$, $O(n^4)$ when $p$ is 2 or 3, $O(n^6)$ when $p$ is 4 or 5, and so on.
n = [floor(Int,10^y) for y in LinRange(1, 2, 10)]
error = zeros(10,7)
f = (x,p) -> x + x.^p.*(2-x).^p
for p ∈ 1:7,
S = trapezoidal(x->f(x,p),[0,2],10^6)
for i ∈ 1:length(n)
Sₙ = trapezoidal(x->f(x,p),[0,2],n[i])
error[i,p] = abs(Sₙ - S)/S
end
end
slope = ([log.(n) one.(n)]\log.(error))[1,:]
info = .*(string.((1:7)'),": slope=",string.(round.(slope')))
plot(n,error, xaxis=:log, yaxis=:log, labels = info)
Clenshaw–Curtis quadrature. applies the trapezoidal rule to a discrete cosine transform (type-1) as a means of numerically evaluating the integral $\int_{-1}^{1} f(x) \,\mathrm{d}x$. We'll test the integral on the function $f(x) = 8 \cos x^2$, with an integral of approximately 0.566566
using FFTW, LinearAlgebra
function clenshaw_curtis(f,n)
x = cos.(π*(0:n)'/n)
w = zeros(n+1,1); w[1:2:n+1] = 2 ./ (1 .- (0:2:n).^2)
1/n * dctI(f.(x)) ⋅ w
end
function dctI(f)
g = FFTW.r2r!([f...],FFTW.REDFT00)
[g[1]/2; g[2:end-1]; g[end]/2]
end
clenshaw_curtis(x -> cos(8x^2),20)
A mathematical comment: we could have also defined a type-1 DCT explicitly in terms of its underlying FFTs if we wanted to crack the black box open just a wee bit more.
function dctI(f)
n = length(f)
g = real(fft(f[[1:n; n-1:-1:2]]))
[g[1]/2; g[2:n-1]; g[n]/2]
end
Gauss–Legendre quadrature. We first compute the Legendre weights and nodes and then apply Gauss–Legendre quadrature to compute $$\int_{-1}^{1} \cos x \cdot \mathrm{e}^{-x^2} \,\mathrm{d}x$$ using a nominal number of nodes $n=8$.
function gauss_legendre(n)
a = zeros(n)
b = (1:n-1).^2 ./ (4*(1:n-1).^2 .- 1)
𝟙² = 2
λ, v = eigen(SymTridiagonal(a, sqrt.(b)))
(λ, 𝟙²*v[1,:].^2)
end
n = 8
f = x -> cos(x)*exp(-x^2);
nodes, weights = gauss_legendre(n)
f.(nodes) ⋅ weights
Alternatively, the Julia library FastGaussQuadrature.jl provides a fast, accurate method of computing the nodes and weights
using FastGaussQuadrature
nodes, weights = gausslegendre(n)
f.(nodes) ⋅ weights
r = exp.(2im*π*(0:.01:1))
plot(@. (1.5r^2 - 2r + 0.5)/r^2); plot!(aspect_ratio=:equal)
Multistep coefficients. The function multistepcoefficients
determines the multistep coefficients for stencil given by m
and n
. The function plotstability
uses these coefficients to plot boundary of the region of absolute stability. We'll test it on the Adams–Moulton method with input m = [0 1]
and n = [0 1 2]
.
function multistepcoefficients(m,n)
s = length(m) + length(n) - 1
A = (m.+1).^(0:s-1) .|> Rational
B = (0:s-1).*((n.+1).^[0;0:s-2])
c = -[A[:,2:end] B]\ones(Int64,s)
[1;c[1:length(m)-1]], c[length(m):end]
end
function plotstability(a,b)
λk(r) = (a ⋅ r.^-(0:length(a)-1)) ./ (b ⋅ r.^-(0:length(b)-1))
r = exp.(im*LinRange(0,2π,200))
plot(λk.(r),label="",aspect_ratio=:equal)
end
m = [0 1]; n = [0 1 2]
a, b = zeros(maximum(m)+1), zeros(maximum(n)+1)
a[m.+1], b[n.+1] = multistepcoefficients(m,n)
plotstability(a,b)
Recipe for solving an ODE. The general recipe for solving an ODE is to
Let's apply this recipe to solve the pendulum problem $u'' = \sin u$ with initial conditions $u(0) = \pi/9$ and $u'(0) = 0$ over $t\in[0,8\pi]$.
using DifferentialEquations, Plots
pendulum(u,p,t) = [u[2]; -sin(u[1])]
u₀ = [8π/9,0]; tspan = [0,8π]
problem = ODEProblem(pendulum, u₀, tspan)
method = Trapezoid()
solution = solve(problem,method)
plot(solution, xaxis="t", label=["θ" "ω"])
ODE benchmarking. Let's benchmark the efficiencies of several numerical ODE solvers on the Lotka–Volterra equations. We'll first generate a high-accuracy solution to approximate the exact solution.
using DifferentialEquations, ParameterizedFunctions, DiffEqDevTools, Sundials
f = @ode_def LotkaVolterra begin
dx = α*x - β*x*y
dy = -γ*y + δ*x*y
end α β γ δ
params = [1.5,1.0,3.0,1.0]
problem = ODEProblem(f,[1.0;1.0],(0.0,10.0),params)
solution = solve(problem,Vern7(),abstol=1e-15,reltol=1e-15)
testsol = TestSolution(solution)
plot(solution)
Now, we'll choose a set of solvers to benchmark.
solvers = [
Dict(:alg=>VCABM())
Dict(:alg=>DP5())
Dict(:alg=>BS3())
Dict(:alg=>CVODE_Adams())
Dict(:alg=>Tsit5())
Dict(:alg=>DP8())
Dict(:alg=>Vern8())
Dict(:alg=>Vern9())
];
We use DiffEqDevTools.jl to generate a working precision set.
abstols = 10 .^[-(LinRange(6,13,40))...]
reltols = 1e3*abstols
wps = WorkPrecisionSet(problem,abstols,reltols,solvers;
appxsol=testsol,save_everystep=true,numruns=100,maxiters=10^5)
Let's plot the data and fit regression lines. First, we'll remove the outliers.
using DataFrames, LinRegOutliers
function logfit(w, c = 1)
df = DataFrame(x=log.(w.errors), y=log.(w.times))
reg = createRegressionSetting(@formula(y ~ x), df)
s = lts(reg; crit = c)
β = s["betas"]
delete!(df, s["outliers"])
x = exp.([minimum(df.x) maximum(df.x)])
y = [exp(β[1])*(x[1])^β[2] exp(β[1])*(x[2])^β[2]]
(x,y,exp.(df.x),exp.(df.y))
end
Now, we plot the compute time versus precision.
plot(xaxis=:log,yaxis=:log,xflip = true)
for i∈1:length(wps)
x, y, xₚ, yₚ = logfit(wps[i],5)
scatter!(xₚ,yₚ,alpha=0.25,labels=:none,markersize=5,color=i)
plot!(x',y',linewidth=3,labels=:none,color=i)
annotate!(x[1], y[1], text.(wps.names[i],9,:left))
end
plot!()
Differential-algebraic equations. The pendulum problem.
using DifferentialEquations
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
U₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0]
function pendulum(dU,U,p,t)
x,y,u,v = U; ℓ,g,m = p
τ = m*(u^2 + v^2 + g*y)/ℓ^2
dU .= (u, v, -τ*x/m, -τ*y/m + g)
end
problem = ODEProblem(pendulum, U₀, tspan, (ℓ,1,1))
solution = solve(problem, Tsit5())
plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)
function residual(r,w,p,t)
x,y,u,v = w; ℓ,g,m = p
r .= (x^2+y^2-ℓ^2, x*u+y*v, 0, 0)
end
cb = ManifoldProjection(residual)
solution = solve(problem, Tsit5(), callback=cb)
plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)
using DifferentialEquations
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, -ℓ*cos(θ₀)]
function pendulum(dU,U,p,t)
x,y,u,v,τ = U; ℓ,g,m = p
dU .= (u, v, -τ*x/m, -τ*y/m + g, -ℓ^2*τ + m*g*y + m*(u^2 + v^2))
end
M = Diagonal([1,1,1,1,0])
f = ODEFunction(pendulum, mass_matrix=M)
problem = ODEProblem(f,u₀,tspan,(ℓ,1,1))
solution = solve(problem, Rodas5(), reltol=1e-8, abstol=1e-8)
plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)
using ModelingToolkit, DifferentialEquations
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, 0]
M = Diagonal([1,1,1,1,0])
function pendulum(dU,U,p,t)
x,y,u,v,τ = U; ℓ,g,m = p
dU .= (u, v, -τ*x/m, -τ*y/m + g, x^2 + y^2 - ℓ^2)
end
f = ODEFunction(pendulum, mass_matrix=M)
problem = ODEProblem(f, u₀, tspan, [ℓ,1,1])
sys = modelingtoolkitize(problem)
pendulum_sys = structural_simplify(dae_index_lowering(sys))
problem = ODAEProblem(pendulum_sys, [], tspan)
solution = solve(problem, Tsit5(), abstol=1e-8, reltol=1e-8)
plot(solution, idxs=(1,3), yflip=true, aspect_ratio=:equal)
Gauss–Legendre–Lobatto orthogonal collocation The following code solves the pendulum problem using orthogonal collocation.
using FastGaussQuadrature
function differentiation_matrix(n,Δt=1)
nodes, _ = gausslobatto(n+1)
t = (nodes[2:end].+1)/2
A = t.^(0:n-1)'.*(1:n)'
B = t.^(1:n)'
(A/B)/Δt, t
end
function pendulumGLL(U,(U₀,D,n))
x,y,u,v,τ = U[1:n],U[n+1:2n],U[2n+1:3n],U[3n+1:4n],U[4n+1:5n]
x₀,y₀,u₀,v₀,_ = U₀
ℓ,g,m = (1,1,1)
[D(x,x₀) .- u;
D(y,y₀) .- v;
D(u,u₀) .+ τ.*x/m;
D(v,v₀) .+ τ.*y/m .- g;
x.^2 + y.^2 .- ℓ^2]
end
using NonlinearSolve
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, 0]
n = 5; N = 100; Δt = 30/N
M,_ = differentiation_matrix(n,Δt)
D = (u,u₀) -> M*(u .- u₀)
u = u₀
for i = 2:N
problem = NonlinearProblem(pendulumGLL,[ones(n)*u₀'...],(u₀,D,n))
solution = solve(problem,NewtonRaphson(),abstol=1e-12)
u = [u reshape(solution.u,5,n)']
u₀ = u[:,end]
end
plot(u[1,:], u[2,:], yflip=true, aspect_ratio=:equal)
Heat equation using the backward Euler method. Let's solve the heat equation using the backward Euler method with initial conditions given by a rectangular function and absorbing boundary conditions.
Δx = .01; Δt = .01; L = 2; ν = Δt/Δx^2; uₗ = 0; uᵣ = 0;
x = -L:Δx:L; m = length(x)
u = (abs.(x).<1)
u[1] += 2ν*uₗ; u[m] += 2ν*uᵣ
D = Tridiagonal(ones(m-1), -2ones(m), ones(m-1))
D[1,2] = 0; D[m,m-1] = 0
A = I - ν*D
for i in 1:20
u = A\u
end
plot(x,u)
Heat equation using the Crank–Nicolson method. Let's solve the heat equation again using the Crank–Nicolson method with initial conditions given by a rectangular function. This time we'll use reflecting boundary conditions. Notice how the high-frequency information does not decay as it did when using the backward Euler method.
Δx = .01; Δt = .01; L = 2; ν = Δt/Δx^2
x = -L:Δx:L; m = length(x)
u = float.(abs.(x).<1)
D = Tridiagonal(ones(m-1), -2ones(m), ones(m-1))
D[1,2] = 2; D[m,m-1] = 2
A = 2I + ν*D
B = 2I - ν*D
anim = @animate for i in 1:40
plot(x,u, legend=:none, ylims=(0,1))
u .= B\(A*u)
end
gif(anim, "heat_cn.gif", fps = 5)
Porous medium equation. We'll now solve the porous medium equation $u_t = (u^2u_x)_x$ using an adaptive-step BDF routine.
using Sundials
m = 400; L = 2; x = LinRange(-L,L,m); Δx = x[2]-x[1]
α = (u -> u.^2)
Du(u,Δx,t) = [0;diff(α((u[1:end-1]+u[2:end])/2).*diff(u))/Δx^2;0]
u₀ = (abs.(x).<1)
problem = ODEProblem(Du,u₀,(0,2),Δx)
method = CVODE_BDF(linear_solver=:Band, jac_lower=1, jac_upper=1)
solution = solve(problem, method);
anim = @animate for t in LinRange(0,2,200)
plot(x,solution(t),legend=:none,fill=(0,0.4,:red),ylims=(0,1))
end
gif(anim, "porous.gif", fps = 15)
using Interact
@manipulate for t∈slider(0:0.01:2; value=0, label="time")
plot(x,solution(t), fill = (0, 0.4, :red))
plot!(ylims=(0,1),legend=:none)
end
Heat equation. The formal solution to the heat equation is $$u(t,x) = \mathrm{F}^{-1}\left[ \mathrm{e}^{-\xi^2 t} \mathrm{F} u(0,x) \right].$$
using FFTW
m = 256; ℓ = 4
ξ² = (fftshift(-(m-1)÷2:m÷2)*(2π/ℓ)).^2
u = (t,u₀) -> real(ifft(exp.(-ξ²*t).*fft(u₀)))
x = LinRange(-ℓ/2,ℓ/2,m+1)[2:end]
u₀ = (abs.(x).<1)
plot(x,u(0.1,u₀))
F = plan_fft(u₀)
u = (t,u₀) -> real(F\(exp.(-ξ²*t).*(F*u₀)))
Kuramoto–Sivashinsky equation. Let's solve the KSE using an IMEX method.
#= In moving from Julia 1.7 to 1.9 this code no longer works.
using FFTW, DifferentialEquations, SciMLOperators
ℓ = 100; T = 200.0; n = 512
x = LinRange(-ℓ/2,ℓ/2,n+1)[2:end]
u₀ = exp.(-x.^2)
iξ = im*[0:(n÷2);(-n÷2+1):-1]*2π/ℓ
F = plan_fft(u₀)
L = DiagonalOperator(-iξ.^2-iξ.^4)
N = (u,_,_) -> -0.5iξ.*(F*((F\u).^2))
problem = SplitODEProblem(L,N,F*u₀,(0,T))
method = KenCarp3(linsolve=KrylovJL(),autodiff=false)
solution = solve(problem,method,reltol=1e-12)
u = t -> real(F\solution(t))
=#
#=
using ImageShow,ColorSchemes
s = hcat(u.(LinRange(0,T,500))...)
get(colorschemes[:binary],abs.(s),:extrema)
=#
Navier–Stokes equation. We solve the Navier–Stokes equation in three parts: define the functions, initialize the variables, and iterate over time.
Δᵒ(Q,step=1) = Q - circshift(Q,(step,0))
flux(Q,c) = c.*Δᵒ(Q) - 0.5c.*(1 .- c).*(Δᵒ(Q,1)+Δᵒ(Q,-1))
H = (u,v,iξ₁,iξ₂) -> F*((F\u).*(F\(iξ₁.*u)) + (F\v).*(F\(iξ₂.*u)))
using FFTW
ℓ = 2; n = 128; ϵ = 0.001; Δt = 0.001; Δx = ℓ/n
x = LinRange(Δx,ℓ,n)'; y = x'
Q = (@. 0.5(1 + tanh(10 - 20abs(1 - 2y/ℓ)))).*ones(1,n)
uⁿ = Q .* (1 .+ 0.5sin.(ℓ*π*x))
vⁿ = zeros(n,n)
F = plan_fft(uⁿ)
uⁿ, vⁿ = F*uⁿ, F*vⁿ
uᵒ, vᵒ = uⁿ, vⁿ
iξ₁ = im*fftfreq(n,n)'*(2π/ℓ)
iξ₂ = transpose(iξ₁)
ξ² = iξ₁.^2 .+ iξ₂.^2
H₁ⁿ, H₂ⁿ = H(uⁿ,vⁿ,iξ₁,iξ₂), H(vⁿ,uⁿ,iξ₂,iξ₁)
M₊ = (1/Δt .+ (ϵ/2)*ξ²)
M₋ = (1/Δt .- (ϵ/2)*ξ²);
for i = 1:1200
Q -= flux(Q,(Δt/Δx).*real(F\vⁿ)) + flux(Q',(Δt/Δx).*real(F\uⁿ)')'
H₁ⁿ⁻¹, H₂ⁿ⁻¹ = H₁ⁿ, H₂ⁿ
H₁ⁿ, H₂ⁿ = H(uⁿ,vⁿ,iξ₁,iξ₂), H(vⁿ,uⁿ,iξ₂,iξ₁)
uᵒ = uⁿ - uᵒ + (-1.5H₁ⁿ + 0.5H₁ⁿ⁻¹ + M₊.*uⁿ)./M₋
vᵒ = vⁿ - vᵒ + (-1.5H₂ⁿ + 0.5H₂ⁿ⁻¹ + M₊.*vⁿ)./M₋
ϕ = (iξ₁.*uᵒ + iξ₂.*vᵒ)./(ξ² .+ (ξ².≈ 0))
uⁿ, vⁿ = uᵒ - iξ₁.*ϕ, vᵒ - iξ₂.*ϕ
end
using ColorSchemes
imresize(get(ColorSchemes.acton, clamp01.(Q)),ratio=3)
1.4. Invertibility of random (0,1) matrices. The number of invertible $n\times n$ (0,1) matrices is known for $n$ up to 8. (See the On-Line Encyclopedia of Integer Sequences.) We'll approximate the ratio of invertible matrices by checking the determinants of randomly drawn ones. Let's also plot the known ratio for $n=1,2,\dots,8$.
using LinearAlgebra, Plots
N = 10000
n = [sum([!(det(rand(Bool,d,d))≈0)
for i in 1:N]) for d in 1:20]
plot(n/N,marker=:o)
p = [1,6,174,22560,12514320,28836612000,270345669985440,10160459763342013440]
r = @. 2^(log2(p)-(1:8)^2)
plot!(r,marker=:square,opacity=0.5)
A = {{0,1,0,1},{1,0,2,1},{0,2,0,0},{1,1,0,0}};
W = Inverse[IdentityMatrix[4] - z*A];
Simplify[W[[1,3]]]
Series[%,{z,0,10}]
2.4. Good Will Hunting problem. The following function computes the number of walks for $n$-step walks for the graph
A = [0 1 0 1;1 0 2 1;0 2 0 0;1 1 0 0]
walks(A,i,j,N) = [(A^n)[i,j] for n in 1:N]
walks(A,1,3,9)
2.3. Naive algorithm for the determinant. The determinant of matrix $\mathbf{A}$ is given by the product of the elements along the diagonal of $\mathbf{U}$ multiplied by the parity of the permutation matrix $\mathbf{P}$ from the LU decomposition $\mathbf{PLU} = \mathbf{A}$.
function detx(A)
L,U,P = lu(A)
s = 1
for i in 1:length(P)
m = findfirst(P.==i)
i!=m && ( s *= -1; P[[m,i]]=P[[i,m]] )
end
s * prod(diag(U))
end
A = rand(20,20)
detx(A) ≈ det(A)
2.5. Reverse Cuthill–McKee algorithm. The following function implements a naive reverse Cuthill–McKee algorithm for symmetric matrices. We'll verify the algorithm by applying it to a sparse, random (0,1) matrix.
function rcuthillmckee(A)
r = sortperm(vec(sum(A.!=0,dims=2)))
p = Int64[]
while ~isempty(r)
q = [popfirst!(r)]
while ~isempty(q)
q₁ = popfirst!(q)
append!(p,q₁)
k = findall(!iszero, A[q₁,r])
append!(q,splice!(r,k))
end
end
reverse(p)
end
using SparseArrays
A = sparse(Symmetric(sprand(1000,1000,0.003)))
p = rcuthillmckee(A)
plot(spy(A), spy(A[p,p]), colorbar = false)
2.6. Dolphins of Doubtful Sound. We'll reuse the code above used to draw the original graph of the dolphins.
using DelimitedFiles, SparseArrays, Graphs, GraphPlot
ij = readdlm(download(bucket*"dolphins.csv"),',',Int)
A = sparse(ij[:,1],ij[:,2],one.(ij[:,1]))
gplot(Graph(A),layout=circular_layout)
p = rcuthillmckee(A)
gplot(Graph(A[p,p]),layout=circular_layout)
2.8. Stigler diet problem. Let's solve the Stigler diet problem. We'll use the naïve simplex
function defined above.
using CSV, DataFrames
diet = DataFrame(CSV.File(download(bucket*"diet.csv")))
A = Array(diet[2:end,4:end])'
c = ones.(size(A,2))
b = Array(diet[1,4:end])
food = diet[2:end,1]
solution = simplex(b,A',c)
print("foods: ", food[solution.y .!= 0], "\n")
print("daily cost: ", solution.z)
using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, x[1:size(A,2)] ≥ 0)
@objective(model, Min, c' * x)
@constraint(model, A * x .≥ b)
optimize!(model)
print("foods: ", food[JuMP.value.(x) .!= 0], "\n")
print("daily cost: ", objective_value(model))
2.9. Six degrees of Kevin Bacon. Let's determine the shortest path between two actors along with their connecting movies. We'll first define helper functions get_names
and get_adjacency_matrix
. Then we'll build an biadjacency matrix $\mathbf{B}$ and construct the adjacency matrix $$\mathbf{A} = \begin{bmatrix} \mathbf{0} & \mathbf{B}^\mathsf{T} \\ \mathbf{B} & \mathbf{0}\end{bmatrix}.$$
function findpath(A,a,b)
r = collect(1:size(A,2))
q = [a]; p = [-1]; i = 0
splice!(r,a)
while i<length(q)
k = findall(!iszero, A[q[i+=1],r])
any(r[k].==b) && return(append!(q[backtrack(p,i)],b))
append!(q,splice!(r,k))
append!(p,fill(i,length(k)))
end
display("No path.")
end
function backtrack(p,i)
s = []; while (i!=-1); append!(s,i); i = p[i]; end
return(reverse(s))
end
using DelimitedFiles, SparseArrays
get_names(filename) = readdlm(download(bucket*filename*".txt"),'\n')
function get_adjacency_matrix(filename)
i = readdlm(download(bucket*filename*".csv"),',',Int,'\n')
sparse(i[:,1], i[:,2], 1)
end
actors = get_names("actors")
movies = get_names("movies")
B = get_adjacency_matrix("actor-movie");
(m,n) = size(B)
A = [spzeros(n,n) B' ; B spzeros(m,m)];
actormovie = [ actors ; movies ];
index = x -> findfirst(actors.==x)[1]
actormovie[findpath(A, index("Emma Watson"), index("Kevin Bacon"))]
actormovie[findpath(A, index("Bruce Lee"), index("Tommy Wiseau"))]
2.10. Sparse matrix storage efficiency.
d = 0.5; A = sprand(800,600,d)
Base.summarysize(A)/Base.summarysize(Matrix(A))
3.4. Image deblurring. Take X
to grayscale image (with pixel intensity between 0 and 1), A
and B
to be Gaussian blurring matrices with standard deviations of 20 pixels, and N
to be a matrix of random values from the uniform distribution over the interval (0,0.01). We'll compare three deblurring methods: inverse, Tikhonov regulation, and the pseuodinverse. We can find a good value for regulariation parameter α = 0.05 with some trial-and-error.
X = load(download(bucket*"laura.png")) .|> Gray
m,n = size(X)
blur = x -> exp(-(x/20).^2/2)
A = [blur(i-j) for i=1:m, j=1:m]; A ./= sum(A,dims=2)
B = [blur(i-j) for i=1:n, j=1:n]; B ./= sum(B,dims=1)
N = 0.01*rand(m,n)
Y = A*X*B + N;
α = 0.05
X₁ = A\Y/B
X₂ = (A'*A+α^2*I)\A'*Y*B'/(B*B'+α^2*I)
X₃ = pinv(A,α)*Y*pinv(B,α)
Gray.([X Y X₁ X₂ X₃])
3.5. Filippelli problem. The National Institute for Standards and Technology (NIST) contrived the Filippelli dataset to benchmark linear regression software. The Filippelli problem consists of fitting a 10th degree polynomial to the data set—a rather ill-conditioned problem. We first need to download the data. Then we'll define three methods for solving the Vandermonde problem: the normal equation, QR decomposition, and the pseudoinverse.
vandermonde(t,n) = vec(t).^(0:n-1)'
build_poly(c,X) = vandermonde(X,length(c))*c
function solve_filip(x,y,n)
V = vandermonde(x, n)
c = Array{Float64}(undef, 3, n)
c[1,:] = (V'*V)\(V'*y)
Q,R = qr(V)
c[2,:] = R\(Matrix(Q)'*y)
c[3,:] = pinv(V,1e-9)*y
r = [norm(V*c[i,:]-y) for i in 1:3]
return(c,r)
end
Now, let's solve the problem and plot the results. Let's also list the coefficients from each method alongside the official NIST coefficients. What do you notice about the coefficients? What methods performs the best?
using DelimitedFiles
data = readdlm(download(bucket*"filip.csv"),',',Float64)
coef = readdlm(download(bucket*"filip-coeffs.csv"),',')
(x,y) = (data[:,2],data[:,1])
β,r = solve_filip(x, y, 11)
X = LinRange(minimum(x),maximum(x),200)
Y = [build_poly(β[i,:],X) for i in 1:3]
plot(X,Y); scatter!(x,y,opacity=0.5)
[coef β']
What makes the Filippelli problem difficult is that the system's huge condition number. We can reduce the condition number by first standardizing the data before using it—i.e., subtracting the mean and dividing by the standard deviation. Look at the difference in condition numbers of the Vandermonde matrix before and after standardizing the data.
using Statistics
zscore(X,x=X) = (X .- mean(x))/std(x)
c,r = solve_filip(zscore(x), zscore(y), 11)
Y = [build_poly(c[i,:],zscore(X,x))*std(y).+mean(y) for i in 1:3]
plot(X,Y); scatter!(x,y,opacity=0.5)
3.7. Modeling daily temperatures. We'll use $u(t) = c_1 \sin(2\pi t) + c_2 \cos(2\pi t) + c_3$ to model the daily temperatures using data recorded in Washington, DC. between 1967 and 1971.
using CSV, DataFrames, Dates
data = download(bucket*"dailytemps.csv")
df = DataFrame(CSV.File(data))
day = Dates.value.(df[:,:date] .- df[1,:date])/365
u = df[:,:temperature]
tempsmodel(t) = [sin.(2π*t) cos.(2π*t) one.(t)]
c = tempsmodel(day)\u
scatter(df[:,:date],u,alpha=0.3)
plot!(df[:,:date],tempsmodel(day)*c,width=3)
3.8. Image recognition. We'll practice identifying handwritten digits using the MNIST database.
using MLDatasets, Arpack, Images
image_train, label_train = MNIST(split=:train)[:]
image_train = reshape(permutedims(image_train,[3,2,1]),60000,:)
V = [( D = image_train[label_train.==i,:];
svds(D,nsv=12)[1].V ) for i ∈ 0:9];
pix = -abs.(reshape(V[3+1],28,:))
rescale = scaleminmax(minimum(pix), maximum(pix))
pix .|> rescale .|> Gray
image_test, label_test = MNIST(split=:test)[:]
image_test = reshape(permutedims(image_test,[3,2,1]),10000,:)
r = [( q = image_test*V[i]*V[i]' .- image_test;
sum(q.^2,dims=2) ) for i∈1:10]
r = hcat(r...)
prediction = [argmin(r[i,:]) for i∈1:10000] .- 1
[sum(prediction[label_test.== i].== j) for j∈0:9, i∈0:9]
3.9. Actor similarity model. We use SVD to find a lower-dimensional subspace relating actors and genres. Then we find the closest actors in that subspace using cosine similarity. We use the functions get_names
and get_adjacency_matrix
developed earlier.
actors = get_names("actors")
genres = get_names("genres")
A = get_adjacency_matrix("movie-genre")
A = (A*Diagonal(1 ./ sum(A,dims=1)[:]))
B = get_adjacency_matrix("actor-movie");
using Arpack
X,_ = svds(A*B,nsv=10)
U = X.U; Σ = Diagonal(X.S); Vᵀ = X.Vt
Q = Vᵀ./sqrt.(sum(Vᵀ.^2,dims=1));
index = x -> findfirst(actors.==x)[1]
q = Q[:,index("Steve Martin")][:]
z = Q'*q
r = sortperm(-z)
actors[r[1:10]]
p = U*Σ*q
r = sortperm(-p)
[genres[r] p[r]/sum(p)]
3.10. Multilateration. We use ordinary least squares and total least squares to solve a multilateration problem. The function tls
is defined above.
X = [3 3 12; 1 15 14; 10 2 13; 12 15 14; 0 11 12]
reference = X[1,:]; X .-= reference'
A = [2 2 -2].*X; b = (X.^2)*[1; 1; -1]
xₒₗₛ = A\b + reference
xₜₗₛ = tls(A,b) + reference
xₒₗₛ, xₜₗₛ
3.11. ShotSpotter. Let's modify the previous solution for this multilateration problem. We'll first download the data set.
using CSV, DataFrames
df = DataFrame(CSV.File(download(bucket*"shotspotter.csv")))
X = Array(df[1:end-1,:]); xₒ = Array(df[end,:]);
c² = 328.9^2
reference = X[1,:]; X .-= reference'
A = [2 2 2 -2c²].*X; b = (X.^2)*[1; 1; 1; -c²]
xₒₗₛ = A\b + reference
error = xₒₗₛ - xₒ
4.1. Girko's circular law. Let's plot out the eigenvalues of a few thousand normal random matrices of size $n$ to get a probability distribution in the complex plane. What do you notice?
n = 8
E = collect(Iterators.flatten([eigvals(randn(n,n)) for i∈1:2000]))
scatter(E,mc=:black,ma=0.05,legend=nothing,aspect_ratio=:equal)
4.4. Rayleigh quotient iteration. Let's define a function that finds an eigenvalue $\lambda$ and eigenvector $\mathbf{x}$ of a matrix. The function will choose a random initial guess for $\mathbf{x}$ unless one is provided. We'll then verify the algorithm on a symmetric matrix. Rayleigh quotient iteration works on general classes of matrices, but it often has difficulty converging when matrices get large or far from symmetric—i.e., when eigenvectors get closer together.
function rayleigh(A)
x = randn(size(A,1),1)
while true
x = x/norm(x)
ρ = (x'*A*x)[1]
M = A - ρ*I
abs(cond(M, 1)) < 1e12 ? x = M\x : return(ρ,x)
end
end
A = [2 3 6 4; 3 0 3 1; 6 3 8 8; 4 1 8 2]
rayleigh(A)
4.5. Implicit QR method. We'll define a function that computes all the eigenvalues of a matrix using the implicit QR method. We'll then verify the algorithm on a matrix with known eigenvalues.
function implicitqr(A)
n = size(A,1)
tolerance = 1E-12
H = Matrix(hessenberg(A).H)
while true
if abs(H[n,n-1])<tolerance
if (n-=1)<2; return(diag(H)); end
end
Q = givens([H[1,1]-H[n,n];H[2,1]],1,2)[1]
H[1:2,1:n] = Q*H[1:2,1:n]
H[1:n,1:2] = H[1:n,1:2]*Q'
for i = 2:n-1
Q = givens([H[i,i-1];H[i+1,i-1]],1,2)[1]
H[i:i+1,1:n] = Q*H[i:i+1,1:n]
H[1:n,i:i+1] = H[1:n,i:i+1]*Q'
end
end
end
n = 20; S = randn(n,n)
D = Diagonal(1:n); A = S*D*inv(S)
implicitqr(A)
4.6. Hollywood eigenvector centrality. We'll use eigenvector centrality to determine who's at the center of Hollywood. We use the functions get_names
and get_adjacency_matrix
developed earlier.
actors = get_names("actors")
B = get_adjacency_matrix("actor-movie")
M = (B'*B .!= 0) - I
v = ones(size(M,1))
for k = 1:8
v = M*v; v /= norm(v);
end
r = sortperm(-v); actors[r][1:10]
What is Kevin Bacon's ranking?
findfirst(actors[r] .== "Kevin Bacon")
Let's also compute degree centrality, which ranks each node by its number of edges, i.e., the number of actors with which an actor has appeared.
actors[ -sum(M,dims=1)[:] |> sortperm ][1:10]
4.7. Randomized SVD. We define a method that implements randomized SVD. The idea is to start with a set of $k$ random vectors and perform a few steps of the naive QR method to generate a $k$-dimensional subspace closer to the space of dominant singular values. Then, we perform SVD on that subspace. We may not get the exact singular values, but we just need a good guess. Because the matrix is huge, this method can be significantly faster than SVD or sparse SVD.
function randomizedsvd(A,k)
Z = rand(size(A,2),k);
Q = Matrix(qr(A*Z).Q)
for i in 1:3
Q = Matrix(qr(A'*Q).Q)
Q = Matrix(qr(A*Q).Q)
end
W = svd(Q'*A)
return((Q*W.U,W.S,W.V))
end
using Images
img = load(download(bucket*"red-fox.jpg"))
A = Float64.(Gray.(img))
U,S,V = randomizedsvd(A,10)
Gray.([A U*Diagonal(S)*V'])
using Arpack
@time randomizedsvd(A,10)
@time svds(A, nsv=10)
@time svd(A);
ϵ = 1 .- sqrt.(cumsum(S.^2))/norm(A)
plot(ϵ,marker=:circle)
4.8. 100-dollar, 100-digit challenge. "The infinite matrix $\mathbf{A}$ with entries a₁₁=1, a₁₂ = 1/2, a₂₁ = 1/3, a₁₃ = 1/4, a₂₂ = 1/5, a₃₁ = 1/6, and so on, is a bounded operator on $\ell^2$. What is $\|\mathbf{A}\|_2$?"
using Arpack
m = 5000
A = [1/(1 + (i+j-1)*(i+j)/2 - j) for i∈1:m, j∈1:m]
svds(A,nsv=1)[1].S[1]
5.4. 3D Poisson equation. Let's compare the Jacobi, Gauss–Seidel, SOR, and conjugate gradient methods on the Poisson equation over the unit cube.
using SparseArrays, LinearAlgebra
⊗(x,y) = kron(x,y); ϕ = x -> x-x^2
n = 50 ; x = (1:n)/(n+1); Δx = 1/(n+1)
J = sparse(I, n, n)
D = spdiagm(-1 => ones(n-1), 0 => -2ones(n), 1 => ones(n-1) )
A = ( D⊗J⊗J + J⊗D⊗J + J⊗J⊗D ) / Δx^2
f = [ϕ(x)*ϕ(y) + ϕ(x)*ϕ(z) + ϕ(y)*ϕ(z) for x∈x, y∈x, z∈x][:]
uₑ = [ϕ(x)*ϕ(y)*ϕ(z)/2 for x∈x, y∈x, z∈x][:];
function stationary(A,b,ω=0,n=400)
ϵ = []; u = zero.(b)
P = (1-ω)*sparse(Diagonal(A)) + ω*sparse(LowerTriangular(A))
for i=1:n
u += P\(b-A*u)
append!(ϵ,norm(u - uₑ,1))
end
return(ϵ)
end
function conjugategradient(A,b,n=400)
ϵ = []; u = zero.(b)
p = r = b-A*u
for i=1:n
Ap = A*p
α = (r'*p)/(Ap'*p)
u += α.*p; r -= α.*Ap
β = (r'*Ap)/(Ap'*p)
p = r - β.*p
append!(ϵ,norm(u - uₑ,1))
end
return(ϵ)
end
ϵ = zeros(400,4)
@time ϵ[:,1] = stationary(A,-f,0)
@time ϵ[:,2] = stationary(A,-f,1)
@time ϵ[:,3] = stationary(A,-f,1.9)
@time ϵ[:,4] = conjugategradient(A,-f)
method = ["Jacobi" "Gauss-Seidel" "SOR" "Conjugate Gradient"]
plot(ϵ,yaxis=:log,labels=method,bglegend=RGBA(1,1,1,0.7))
k = 1:120; ([one.(k) k]\log10.(ϵ[k,:]))[2,:]
5.5. 100-dollar, 100-digit challenge. "Let $\mathbf{A}$ be the 20000$\times$20000 matrix whose entries are zero everywhere except for the primes 2, 3, 5, 7,..., 224737 along the main diagonal and the number 1 in all the positions $a_{ij}$ with |i-j|=1,2,4,8,...,16384. What is the (1,1)-entry of $\mathbf{A}^{-1}$?"
using Primes, LinearAlgebra, SparseArrays, IterativeSolvers
n = 20000
d = 2 .^ (0:14); d = [-d;d]
P = Diagonal(primes(224737))
B = [ones(n - abs(d)) for d∈d]
A = sparse(P) + spdiagm(n ,n, (d .=> B)...)
b = zeros(n); b[1] = 1
cg(A, b; Pl=P)[1]
6.1. Radix-3 FFT. The radix-3 FFT is similar to the radix-2 FFT. We'll verify that the code is correct by comparing it with a built-in FFT
function fftx3(c)
n = length(c)
ω = exp(-2im*π/n)
if mod(n,3) == 0
k = collect(0:n/3-1)
u = [transpose(fftx3(c[1:3:n-2]));
transpose((ω.^k).*fftx3(c[2:3:n-1]));
transpose((ω.^2k).*fftx3(c[3:3:n]))]
F = exp(-2im*π/3).^([0;1;2]*[0;1;2]')
return(reshape(transpose(F*u),:,1))
else
F = ω.^(collect(0:n-1)*collect(0:n-1)')
return(F*c)
end
end
using FFTW
v = rand(24,1)
[fft(v) fftx3(v)]
6.2. Fast multiplication. The following function uses FFTs to multiply two large integers (inputted as strings). We'll verify that the algorithm works by multiplying the RSA-129 factors.
using FFTW
function multiply(p_,q_)
p = [parse.(Int,split(reverse(p_),""));zeros(length(q_),1)]
q = [parse.(Int,split(reverse(q_),""));zeros(length(p_),1)]
pq = Int.(round.(real.(ifft(fft(p).*fft(q)))))
carry = pq .÷ 10
while any(carry.>0)
pq -= carry*10 - [0;carry[1:end-1]]
carry = pq .÷ 10
end
n = findlast(x->x>0, pq)
return(reverse(join(pq[1:n[1]])))
end
p = "32769132993266709549961988190834461413177642967992942539798288533"
q = "3490529510847650949147849619903898133417764638493387843990820577"
multiply(p,q), parse(BigInt, p)*parse(BigInt, q)
using Random
p = randstring('0':'9', 100000)
q = randstring('0':'9', 100000)
@time multiply(p,q)
@time parse(BigInt, p)*parse(BigInt, q);
6.3. Fast discrete cosine transform.
function dctII(f)
n = size(f,1)
ω = exp.(-0.5im*π*(0:n-1)/n)
return real(ω.*fft(f[[1:2:n; n-mod(n,2):-2:2],:],1))
end
function idctII(f)
n = size(f,1)
ω = exp.(-0.5im*π*(0:n-1)/n)
f[1,:] = f[1,:]/2
f[[1:2:n; n-mod(n,2):-2:2],:] = 2*real(ifft(f./ω,1))
return f
end
dct2(f) = dctII(dctII(f')')
idct2(f) = idctII(idctII(f')')
6.4. DCT image compression.
function dctcompress2(A,d)
B = dct(Float64.(A))
m,n = size(A)
m₀,n₀ = sqrt(d).*(m,n) .|> floor .|> Int
B₀ = B[1:m₀,1:n₀]
A₀ = idct([B₀ zeros(m₀,n-n₀); zeros(m-m₀,n)])
A₀ = A₀ |> clamp01! .|> Gray
sizeof(B₀)/sizeof(B), sqrt(1-(norm(B₀)/norm(B))^2), A₀
end
A = Gray.(load(download(bucket*"laura_square.png")))
[A dctcompress2(A,0.05)[3]]
8.7. Aitken $\Delta^2$ process. Let's approximate the Leibniz series $$1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \dots = \pi$$ with partial sums and with Aitken's extrapolation of the partial sums. We'll plot error to examine convergence.
aitken₁(x₁,x₂,x₃) = x₃ - (x₃-x₂)^2/(x₃ - 2x₂ + x₁)
aitken₂(x₁,x₂,x₃) = (x₁*x₃ - x₂^2)/(x₃ - 2x₂ + x₁)
n = 60000
p = cumsum([(-1)^i*4/(2i+1) for i=0:n])
p₁ = aitken₁.(p[1:n-2],p[2:n-1],p[3:n])
p₂ = aitken₂.(p[1:n-2],p[2:n-1],p[3:n])
plot(abs.(π.-p)); plot!(abs.(π.-p₂)); plot!(abs.(π.-p₁))
plot!(xaxis=:log,yaxis=:log)
8.12. Solving a nonlinear system. We'll solve $$(x^2+y^2)^2 - 2 (x^2 - y^2) =0$$ $$(x^2+y^2 -1)^3-x^2y^3 = 0$$ using homotopy continuation and Newton's method.
using DifferentialEquations
function homotopy(f,df,x)
dxdt(x,p,t) = -df(x)\p
sol = solve(ODEProblem(dxdt,x,(0.0,1.0),f(x)))
sol.u[end]
end
function newton(f,df,x)
for i in 1:100
Δx = -df(x)\f(x)
norm(Δx) > 1e-8 ? x += Δx : return(x)
end
end
f = z -> ((x,y)=tuple(z...);
[(x^2+y^2)^2-2(x^2-y^2); (x^2+y^2-1)^3-x^2*y^3])
df = z -> ((x,y)=tuple(z...);
[4x*(x^2+y^2-1) 4y*(x^2+y^2+1);
6x*(x^2+y^2-1)^2-2x*y^3 6y*(x^2+y^2-1)^2-3x^2*y^2])
display(homotopy(f,df,[1,1]))
display(newton(f,df,[1,1]))
8.13. Secp256k1 Elliptic curve Diffie–Hellman. We'll first write a function that implements the point addition and point doubling. Then, we implement the double-and-add algorithm. Finally, we test the algorithm using Alice and Bob.
function ⊕(P,Q)
a = 0
r = BigInt(2)^256 - 2^32 - 977
if P[1]==Q[1]
d = invmod(2*P[2], r)
λ = mod((3*powermod(P[1],2,r)+a)*d, r)
else
d = invmod(Q[1]-P[1], r)
λ = mod((Q[2]-P[2])*d, r)
end
x = mod(powermod(λ, 2, r) - P[1] - Q[1], r)
y = mod(-λ*(x-P[1])-P[2], r)
[x;y]
end
Base.isodd(m::BigInt) = ((m&1)==1)
function dbl_add(m,P)
if m > 1
Q = dbl_add(m>>1,P)
return isodd(m) ? (Q⊕Q)⊕P : Q⊕Q
else
return P
end
end
P₁ = big"0x79BE667EF9DCBBAC55A06295CE87
0B07029BFCDB2DCE28D959F2815B16F81798"
P₂ = big"0x483ADA7726A3C4655DA4FBFC0E11
08A8FD17B448A68554199C47D08FFB10D4B8"
P = [P₁; P₂]
m, n = 1829442, 3727472
mP = dbl_add(m,P)
nmP = dbl_add(n,mP)
nP = dbl_add(n,P)
print("Alice's private key:\n$m\n\n")
print("Bob's private key:\n$n\n\n")
print("Alice's public key:\n$mP\n\n")
print("Bob's public key:\n$nP\n\n")
print("Alice and Bob's shared secret key:\n$nmP")
9.2. Periodic parametric splines. We modify the code spline_natural
(above) to make a generate a spine with periodic boundary conditions. The function evaluate_spline
is duplicated from the code above.
function spline_periodic(x,y)
h = diff(x)
d = 6*diff(diff([y[end-1];y])./[h[end];h])
α = h[1:end-1]
β = 2*(h+circshift(h,1))
C = Matrix(SymTridiagonal(β,α))
C[1,end]=h[end]; C[end,1]=h[end]
m = C\d
return([m;m[1]])
end
function evaluate_spline(xᵢ,yᵢ,m,n)
h = diff(xᵢ)
B = yᵢ[1:end-1] .- m[1:end-1].*h.^2/6
A = diff(yᵢ)./h - h./6 .*diff(m)
x = range(minimum(xᵢ),maximum(xᵢ),length=n+1)
i = sum(xᵢ' .≤ x,dims=2)
i[end] = length(xᵢ)-1
y = @. (m[i]*(xᵢ[i+1]-x)^3 + m[i+1]*(x-xᵢ[i])^3)/(6h[i]) +
A[i]*(x-xᵢ[i]) + B[i]
return(x,y)
end
n = 5; nx = 30
x = rand(n); y = rand(n)
x = [x;x[1]]; y = [y;y[1]]
t = [0;cumsum(sqrt.(diff(x).^2 + diff(y).^2))]
X = evaluate_spline(t,x,spline_periodic(t,x),nx*n)
Y = evaluate_spline(t,y,spline_periodic(t,y),nx*n)
scatter(x,y); plot!(X[2],Y[2],legend=false)
9.3. Radial basis functions. Let's examine how a polynomial $y(x) = \sum_{i=0}^n c_i x^i$ compares with Gaussian and cubic radial basis functions $y(x) = \sum_{i=0}^n c_i \phi(x-x_i),$ taking $\phi(x)= \exp(-20x^2)$ and $\phi(x) = |x|^3$ an interpolant of the Heaviside function.
n = 20; N = 200
x = LinRange(-1,1,n)
y = float(x .> 0);
ϕ₁(x,a) = abs.(x.-a).^3
ϕ₂(x,a) = exp.(-20(x.-a).^2)
ϕ₃(x,a) = x.^a
X = LinRange(-1,1,N)
interp(ϕ,a) = ϕ(X,a')*(ϕ(x,a')\y)
Y₁ = interp(ϕ₁,x)
Y₂ = interp(ϕ₂,x)
Y₃ = interp(ϕ₃,(0:n-1))
scatter(x,y,seriestype = :scatter, marker = :o, legend = :none)
plot!(X,[Y₁,Y₂,Y₃],ylim=[-0.5,1.5])
9.4. Collocation. We'll use collocation to solve Bessel's equation. We first define a function to solve general linear boundary value problems. And, then we define a function to interpolate between collocation points.
function collocation_solve(L,f,bc,x)
h = x[2]-x[1]
S = L(x)*([1 -1/2 1/6; -2 0 2/3; 1 1/2 1/6]./[h^2 h 1])'
d = [bc[1]; f(x); bc[2]]
A = Matrix(Tridiagonal([S[:,1];0], [0;S[:,2];0], [0;S[:,3]]))
A[1,1:3], A[end,end-2:end] = [1 4 1]/6, [1 4 1]/6
return(A\d)
end
function collocation_build(c,x,N)
X = LinRange(x[1],x[end],N)
h = x[2] - x[1]
i = Int32.(X .÷ h .+ 1); i[N] = i[N-1]
C = [c[i] c[i.+1] c[i.+2] c[i.+3]]'
B = (x->[(1-x).^3;4-3*(2-x).*x.^2;4-3*(1+x).*(1-x).^2;x.^3]/6)
Y = sum(C.*hcat(B.((X.-x[i])/h)...),dims=1)
return(Array(X),reshape(Y,(:,1)))
end
Now, we can solve the Bessel equation $xu''+u'+xu =0$ with boundary conditions $u(0)=1$ and $u(b)=0$.
using Roots, SpecialFunctions
n = 20; N = 141
L = (x -> [x one.(x) x])
f = (x -> zero.(x) )
b = fzero(besselj0, 11)
x = range(0,b,length=n)
c = collocation_solve(L,f,[1,0],x)
X,Y = collocation_build(c,x,70)
plot(X,[Y besselj0.(X)])
Finally, we'll explore the error and convergence rate.
N = 10*2 .^(1:7); ϵ = []
for n in N
x = LinRange(0,b,n)
c = collocation_solve(L,f,[1,0],x)
X,Y = collocation_build(c,x,n)
append!(ϵ, norm(Y-besselj0.(X)) / n)
end
slope = ([log.(N) one.(N)]\log.(ϵ))[1]
s = "slope = "*string.(round.(slope,digits=2))
plot(N, ϵ, xaxis=:log, yaxis=:log, marker=:o, label=s)
10.3. Fractional derivatives. We'll plot the fractional derivatives for a function.
using FFTW
n = 256; ℓ = 2
x = (0:n-1)/n*ℓ .- ℓ/2
ξ = [0:(n/2-1); (-n/2):-1]*(2π/ℓ)
f₁ = exp.(-16*x.^2)
f₂ = sin.(π*x)
f₃ = x.*(1 .- abs.(x))
deriv(f,p) = real(ifft((im*ξ).^p.*fft(f)))
using Interact
func = Dict("Gaussian"=>f₁,"sine"=>f₂,"spline"=>f₃)
@manipulate for f in togglebuttons(func; label="function"),
p in slider(0:0.01:2; value=0, label="derivative")
plot(x,deriv(f,p),legend=:none)
end
10.4. Handwriting classification. We'll use Flux to train a convolutional neural net using MNIST data and then test the model.
using MLDatasets, Flux
model = Chain(
Conv((5,5), 1=>6, tanh, pad=SamePad()),
MeanPool((2,2)),
Conv((5,5), 6=>16, tanh),
MeanPool((2,2)),
Conv((5,5), 16=>120, tanh),
Flux.flatten,
Dense(120, 84, tanh),
Dense(84, 10))
image_train, label_train = MLDatasets.MNIST(split=:train)[:]
image_test, label_test = MLDatasets.MNIST(split=:test)[:]
image_train = Flux.unsqueeze(image_train, 3)
image_test = Flux.unsqueeze(image_test, 3)
label_train = Flux.onehotbatch(label_train, 0:9)
label_test = Flux.onehotbatch(label_test, 0:9);
loss(x,y) = Flux.Losses.logitcrossentropy(model(x),y)
parameters = Flux.params(model)
data = Flux.Data.DataLoader((image_train,label_train),batchsize=100)
optimizer = ADAM()
using ProgressMeter
@showprogress for epochs = 1:5
Flux.train!(loss, parameters, data, optimizer)
end
accuracy(x,y) = sum(Flux.onecold(x) .== Flux.onecold(y))/size(y,2)
accuracy(model(image_test),label_test)
10.5. Multilateration. Let's modify the previous solution for this multilateration problem. We'll first download the data set.
using LsqFit
ϵ = (x,p) -> @. √((x[:,1]-p[1])^2 + (x[:,2]-p[2])^2) - (x[:,3]-p[3])
x = [3 3 12; 1 15 14; 10 2 13; 12 15 14; 0 11 12]
curve_fit(ϵ,x,zeros(5),zeros(3)).param
using CSV, DataFrames
df = DataFrame(CSV.File(download(bucket*"shotspotter.csv")))
X = Array(df[1:end-1,:]); xₒ = Array(df[end,:]);
using LsqFit
c = 328.6
ϵ = (x,p)->sqrt.(sum([(x[:,i].-p[i]).^2 for i∈1:3]))-c*(x[:,4].-p[4])
xₙₗₛ = curve_fit(ϵ,X,zero(X[:,1]),X[1,:]).param
error = xₙₗₛ - xₒ
Let's also plot the solution. We'll first define a function to plot circles.
function plot_circle(x,y,r)
t = LinRange(0,2π,100);
plot!(x.+r*cos.(t), y.+r*sin.(t),color="black",opacity=0.5);
end
r = 328.6*(X[:,end] .- xₙₗₛ[end])
scatter(X[:,1], X[:,2], color = "black")
[plot_circle(X[i,1],X[i,2],r[i]) for i in 1:length(r)]
scatter!([xₙₗₛ[1]], [xₙₗₛ[2]], color = "red")
plot!(legend=:none,aspect_ratio=:equal)
11.1. Finite difference approximation. Let's find coefficients to the third-order approximation of $f'(x)$ for nodes at $x$, $x+h$, $x+2h$ and $x+3h$. The second row of the solution will tell us the coefficients (first four columns) and the coefficient of truncation (last column).
d = [0,1,2,3]; n = length(d)
C = inv( d.^(0:n-1)' .// factorial.(0:n-1)' )
[C C*d.^n/factorial(n)]
11.2. Richardson extrapolation. The following code is an iterative version of the recursive richardson
function above:
function richardson(f,x,m)
D = []
for i in 1:m
append!(D, ϕ(f,x,2^i))
for j in i-1:-1:1
D[j] = (4^(i-j)*D[j+1] - D[j])/(4^(i-j) - 1)
end
end
D[1]
end
ϕ = (f,x,n) -> (f(x+1/n) - f(x-1/n))/(2/n)
richardson(x->sin(x),0,4)
11.3. Automatic differentiation. Let's extend the Dual class
above by adding methods for division, cosine, and square root to the class definition. We'll also add a few more help functions. I've copied the cells from earlier to the one below.
struct Dual
value
deriv
end
Dual(x) = Dual(x,1)
value(x) = x isa Dual ? x.value : x
deriv(x) = x isa Dual ? x.deriv : 0
Base.:+(u, v) = Dual(value(u)+value(v), deriv(u)+deriv(v))
Base.:-(u, v) = Dual(value(u)-value(v), deriv(u)-deriv(v))
Base.:*(u, v) = Dual(value(u)*value(v),
value(u)*deriv(v)+value(v)*deriv(u))
Base.:sin(u) = Dual(sin(value(u)), cos(value(u))*deriv(u))
Base.:sqrt(u) = Dual(sqrt(value(u)), deriv(u) / 2sqrt(value(u)))
Base.:/(u, v) = Dual(value(u)/value(v),
(value(u)*deriv(v)-value(v)*deriv(u))/(value(v)*value(v)))
Base.:cos(u) = Dual(cos(value(u)), -sin(value(u))*deriv(u))
function get_zero(f,x)
ϵ = 1e-12; δ = 1
while abs(δ) > ϵ
fx = f(Dual(x))
δ = value(fx)/deriv(fx)
x -= δ
end
return(x)
end
Now, we can define Newton's method using this new Dual class and use it to find the zero of $4\sin x + \sqrt{x}$.
get_zero(x->4sin(x)+sqrt(x),4)
We can find a minimum or maximum of $4\sin x + \sqrt{x}$ by modifying Newton's method.
function get_extremum(f,x)
ϵ = 1e-12; δ = 1
while abs(δ)>ϵ
fx = f(Dual(Dual(x)))
δ = deriv(value(fx))/deriv(deriv(fx))
x -= δ
end
return(x)
end
get_extremum(x->4sin(x)+sqrt(x),4)
11.4. Cauchy differentiation formula. Let's compute the sixth derivative of $\mathrm{e}^x(\cos^3 x + \sin^3 x)^{-1}$ evaluated at $x = 0$ using the Cauchy differentiation formula: $$f^{(p)}(a) = \frac{p!}{2\pi\mathrm{i}} \oint_\gamma \frac{f(z)}{(z-a)^{p+1}} \,\mathrm{d}{z}.$$
function cauchyderivative(f, a, p, n = 20, ϵ = 0.1)
ω = exp.(2π*im*(0:(n-1))/n)
factorial(p)/(n*ϵ^p)*sum(@. f(a+ϵ*ω)/ω^p)
end
f = x -> exp(x)/(cos(x)^3 + sin(x)^3)
cauchyderivative(f, 0, 6)
11.5. Gauss–Legendre quadrature. The following function computes the nodes and weights for Gauss–Legendre quadrature by using Newton's method to find the roots of $\mathrm{P_n}(x)$. We'll verify the function on a toy problem.
function gauss_legendre(n)
x = -cos.((4*(1:n).-1)*π/(4n+2))
Δ = one.(x)
dPₙ = 0
while(maximum(abs.(Δ))>1e-16)
Pₙ, Pₙ₋₁ = x, one.(x)
for k ∈ 2:n
Pₙ, Pₙ₋₁ = ((2k - 1)*x.*Pₙ-(k-1)*Pₙ₋₁)/k, Pₙ
end
dPₙ = @. n*(x*Pₙ - Pₙ₋₁)/(x^2-1)
Δ = Pₙ ./ dPₙ
x -= Δ
end
return(x, @. 2/((1-x^2)*dPₙ^2))
end
f = (x-> 2sqrt(1-x^2))
x,w = gauss_legendre(10)
w ⋅ f.(x)
11.7. Fundamental solution to the heat equation. We'll use Gauss–Hermite quadrature to compute the solution to the heat equation $$u(t,x) = \frac{1}{\sqrt{4\pi t}}\int_{-\infty}^\infty u_0(s) \mathrm{e}^{-(x-s)^2/4t} \;\mathrm{d}s.$$
using FastGaussQuadrature
ξ,w = gausshermite(40)
u₀ = x -> sin(x)
u = (t,x) -> w ⋅ u₀.(x.-2sqrt(t)*ξ)/sqrt(π)
x = LinRange(-12,12,100); plot(x,u.(1,x))
11.8. Monte Carlo integration. The following function the volume of an $d$-dimensional sphere using $n$ samples and $m$ trials. We'll use it to verify that error of Monte Carlo integration is $O(1/\sqrt{n})$.
mc_π(n) = sum(sum(rand(2,n).^2,dims=1).<1)/n*4
m = 20; d = []; N = 2 .^ (1:20)
for n ∈ N
append!(d,sum([abs(π - mc_π(n)) for i∈1:m])/m)
end
s = log.(N.^[0 1])\log.(d)
scatter(N,d,xaxis=:log, yaxis=:log)
plot!(N,exp(s[1]).*N.^s[2])
11.10 Orthogonal collocation We'll solve $u'(t) = \alpha u^2$ using the spectral method and pseudospectral method.
using FastGaussQuadrature
function differentiation_matrix(n)
nodes, _ = gausslobatto(n+1)
t = (nodes[2:end].+1)/2
A = t.^(0:n-1)'.*(1:n)'
B = t.^(1:n)'
A/B,t
end
n = 20
M,t = differentiation_matrix(n)
D = (u,u₀) -> M*(u .- u₀)
using NLsolve
u₀ = 1.0; α = 0.9
F = (u,u₀,α) -> D(u,u₀) .- α*u.^2
u = nlsolve(u->F(u,u₀,α),u₀*ones(n)).zero
plot([0;t], [u₀;u], marker=:circle, legend=false)
N = 30; Δt = 1.0/N; n = 8
M,t = differentiation_matrix(n)
D = (u,u₀) -> M*(u .- u₀)/Δt
u₀ = 1.0; U = [u₀]; T = [0.0]; α = 0.9
for i = 0:N-1
u = nlsolve(u->F(u,u₀,α),u₀*ones(n)).zero
u₀ = u[end]
append!(T,(i.+t)*Δt)
append!(U,u)
end
plot(T, U, marker=:circle, legend=false)
12.4. Runge–Kutta stability. The following code plots the region of absolute stability for a Runge–Kutta method with tableau A
and b
.
A = [0 0 0 0 0;
1/3 0 0 0 0;
1/6 1/6 0 0 0;
1/8 0 3/8 0 0;
1/2 0 -3/2 2 0];
b = [1/6 0 0 2/3 1/6];
using LinearAlgebra, Plots
function rk_stability_plot(A,b)
E = ones(length(b),1)
r(λk) = abs.(1 .+ λk * b*((I - λk*A)\E))
x,y = LinRange(-4,4,100),LinRange(-4,4,100)
s = reshape(vcat(r.(x'.+im*y)...),(100,100))
contour(x,y,s,levels=[1],legend=false)
end
rk_stability_plot(A,b)
12.7. Third-order IMEX coefficients. We can determine the coefficients of a third-order IMEX method by inverting a system of linear equations.
i = 0:3
a = ((-i)'.^i.//factorial.(i))\[0;1;0;0]
b = ((-(i.+1))'.^i.//factorial.(i))\[1;0;0;0]
[a b]
12.8. Predictor-corrector stability. We'll use the multistepcoefficients
introduced earlier. The following function provides the orbit of points in the complex plane for an $n$th order Adams–Bashforth–Moulton PE(CE)$^m$.
function multistepcoefficients(m,n)
s = length(m) + length(n) - 1
A = (m.+1).^(0:s-1) .|> Rational
B = (0:s-1).*((n.+1).^[0;0:s-2])
c = -[A[:,2:end] B]\ones(Int64,s)
[1;c[1:length(m)-1]], c[length(m):end]
end
using Polynomials
function PECE(n,m)
_,a = multistepcoefficients([0 1],hcat(1:n-1...))
_,b = multistepcoefficients([0 1],hcat(0:n-1...))
α(r) = a ⋅ r.^(1:n-1)/b[1]
β(r) = b[2:end] ⋅ r.^(1:n-1)/b[1]
z = [(c = [r-1; repeat([r + β(r)],m); α(r)];
Polynomials.roots(Polynomial(c)))
for r in exp.(im*LinRange(0,2π,200))]
hcat(z/b[1]...)
end
z = vcat([PECE(2,i)[:] for i in 0:4]...)
scatter(z,label="",markersize=.5,aspect_ratio=:equal)
12.9. Padé approximant. We'll build a function to compute the coefficients of the Padé approximant to $log(r)$.
function pade(a,m,n)
A = Rational(1)*Matrix(I(m+n+1))
A[:,m+2:end] = [i>j ? -a[i-j] : 0 for i∈1:m+n+1, j∈1:n]
pq = A\a[1:m+n+1]
pq[1:m+1], [1; pq[m+2:end]]
end
m = 3; n = 2
a = [0; ((-1).^(0:m+n).//(1:m+n+1))]
(p,q) = pade(a,m,n)
S = n -> [(-1).^(i+j)*binomial(j,i) for i∈0:n, j∈0:n]
S(m)*p, S(n)*q
12.13. SIR solution. Let's solve the susceptible-infected-recovered (SIR) model for infectious diseases using a general ODE solver.
function SIR!(du,u,p,t)
S,I,R = u
β,γ = p
du[1] = dS = -β*S*I
du[2] = dI = +β*S*I - γ*I
du[3] = dR = +γ*I
end
using DifferentialEquations
u₀ = [0.99; 0.01; 0]
tspan = (0.0,15.0)
p = (2.0,0.4)
problem = ODEProblem(SIR!,u₀,tspan,p)
solution = solve(problem)
plot(solution,labels=["susceptible" "infected" "recovered"])
12.14. Duffing equation. We'll use a high-order, explicit ODE solver for the Duffing equation.
using DifferentialEquations, Plots
function duffing!(dx,x,γ,t)
dx[1] = x[2]
dx[2] = -γ*x[2]+x[1]-x[1]^3+0.3*cos(t)
end
problem = ODEProblem(duffing!,[1.0,0.0],(0.0,200.0),0.37)
solution = solve(problem,Vern7())
plot(solution, idxs=(1,2))
12.15. Shooting method. We'll solve the Airy equation $y'' - xy = 0$ using the shooting method that incorporates an initial value solver into a nonlinear root finder.
function solveIVP(f,u0,tspan)
sol = solve(ODEProblem(f,u0,tspan))
return(sol.u[end][1])
end
using DifferentialEquations, Roots
airy(y,p,x) = [y[2];x*y[1]]
domain = (-12.0,0.0); bc = (1.0,1.0); guess = 5
shoot_airy = (guess -> solveIVP(airy,[bc[1];guess],domain)-bc[2])
v = find_zero(shoot_airy, guess)
Once we have our second initial value, we can plot the solution:
sol = solve(ODEProblem(airy,[bc[1],v],domain))
plot(sol)
13.4. An unconditionally unstable method. Let's generate the coefficients for multistep scheme given by the stencil: We'll use the function multistepcoefficients
introduced earlier. Then we'll plot $r(\lambda k)$ along the real axis. The method is unstable wherever $|r|>1$.
m, n = [0 1 2 3 4], [1]
a, b = zeros(maximum(m)+1), zeros(maximum(n)+1)
a[m.+1], b[n.+1] = multistepcoefficients(m,n)
λk = r -> (a ⋅ r.^-(1:length(a))) ./ (b ⋅ r.^-(1:length(b)))
r = LinRange(0.2,6,100)
plot([λk.(r) λk.(-r)],[r r],xlim=[-2,2])
13.5. Dufort–Frankel method. We'll use the Dufort–Frankel method to solve the heat equation. While this method is unconditionally stable, it generates the wrong solution. Notice that while the long-term behavior is dissipative, the solution is largely oscillatory, and the dynamics are more characteristic of a viscous fluid than heat propagation.
Δx = 0.01; Δt = 0.01
ℓ = 1; x = -ℓ:Δx:ℓ; m = length(x)
Uⁿ = exp.(-8*x.^2); Uⁿ⁻¹ = Uⁿ
ν = Δt/Δx^2; α = 0.5 + ν; γ = 0.5 - ν
B = ν*Tridiagonal(ones(m-1),zeros(m),ones(m-1))
B[1,2] *=2; B[end,end-1] *=2
@gif for i = 1:300
global Uⁿ, Uⁿ⁻¹ = (B*Uⁿ+γ*Uⁿ⁻¹)/α, Uⁿ
plot(x,Uⁿ,ylim=(0,1),label=:none,fill=(0, 0.3, :red))
end
13.7. Schrödinger equation. Let's solve the Schrödinger equation for a harmonic potential using the Strang splitting Crank–Nicolson and confirm that the method is $O(h^2,k^2)$.
function schroedinger(m,n,ε)
x = LinRange(-4,4,m); Δx = x[2]-x[1]; Δt = 2π/n; V = x.^2/2
ψ = exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)
diags = 0.5im*ε*[1 -2 1]/Δx^2 .- im/ε*[0 1 0].*V
D = Tridiagonal(diags[2:end,1], diags[:,2], diags[1:end-1,3])
D[1,2] *= 2; D[end,end-1] *= 2
A = I + 0.5Δt*D
B = I - 0.5Δt*D
for i ∈ 1:n
ψ = B\(A*ψ)
end
return ψ
end
We'll loop over several values for time steps and mesh sizes and plot the error. This may take a while. Go get a snack.
ε = 0.3; m = 20000; eₓ=[]; eₜ=[]
N = floor.(Int,exp10.(LinRange(2,3.7,6)))
x = LinRange(-4,4,m)
ψₘ = -exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)
for n ∈ N
x = LinRange(-4,4,n)
ψₙ = -exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)
append!(eₜ,norm(ψₘ - schroedinger(m,n,ε))/m)
append!(eₓ,norm(ψₙ - schroedinger(n,m,ε))/n)
end
plot(2π./N,eₜ,shape=:circle, xaxis=:log, yaxis=:log)
plot!(8 ./N,eₓ,shape=:circle)
13.8. Polar heat equation. We'll solve a radially symmetric heat equation. Although we divide by zero at $r=0$ when constructing the Laplacian operator, we subsequently overwrite the resulting inf
term when we apply the boundary condition.
t = 0.5; n=100; m = 100
r = LinRange(0,2,m); Δr = r[2]-r[1]; Δt = t/n;
u₀ = @. tanh(32(1 - r)); u = u₀
d = @. [1 -2 1]/Δr^2 + (1/r)*[-1 0 1]/2Δr
D = Tridiagonal(d[2:end,1],d[:,2],d[1:end-1,3])
D[1,1:2] = [-4 4]/Δr^2; D[end,end-1:end] = [2 -2]/Δr^2
A = I - 0.5Δt*D
B = I + 0.5Δt*D
for i = 1:n
u = A\(B*u)
end
plot(r, u,label=:none,fill=(-1, 0.3, :red))
The following is a slower but high-order alternative to the Crank–Nicolson method:
using Sundials
problem = ODEProblem((u,p,t)->D*u,u₀,(0,t))
method = CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)
solution = solve(problem,method)
plot(r, solution(t),label=:none,fill=(-1, 0.3, :red))
13.9. Open boundaries. We can approximate open boundaries by spacing the grid points using a sigmoid function such as $\mathrm{arctanh}\, x$. We start by defining a function logitspace
, a logit analog to np.linspace
. Then we define a Laplacian operator using arbitrary grid spacing. Finally, we solve the heat equation using the Crank–Nicolson using both equally-spaced and logit-spaced grid points.
logitspace(x,n,p) = x*atanh.(LinRange(-p,p,n))/atanh(p)
function laplacian(x)
Δx = diff(x); Δx₁ = Δx[1:end-1]; Δx₂ = Δx[2:end]
d₋ = @. 2/[Δx₁*(Δx₁+Δx₂); Δx₁[1]^2]
d₀ = @. -2/[Δx₂[end].^2; Δx₁*Δx₂; Δx₁[1].^2]
d₊ = @. 2/[Δx₂[end].^2; Δx₂*(Δx₁+Δx₂) ]
Tridiagonal(d₋,d₀,d₊)
end
function heat_equation(x,t,n,u)
Δt = t/n
D² = laplacian(x)
A = I - 0.5Δt*D²
B = I + 0.5Δt*D²
for i ∈ 1:n
u = A\(B*u)
end
return u
end
ϕ = (x,t,s) -> exp.(-s*x.^2/(1+4*s*t))/sqrt(1+4*s*t)
m = 40; t = 15
x = LinRange(-20,20,m)
plot(x,ϕ(x,t,10),label="exact",width=3)
u₁ = heat_equation(x,t,m,ϕ(x,0,10))
plot!(x,u₁,shape=:circle,label="equal")
x = logitspace(20,m,0.999)
u₂ = heat_equation(x,t,m,ϕ(x,0,10))
plot!(x,u₂,shape=:circle,label="logit")
13.10. Allen–Cahn equation. We'll solve the Allen–Cahn equation using Strang splitting.
L = 16; m = 400; Δx = L/m
T = 4; n = 1600; Δt = T/n
x = LinRange(-L/2,L/2,m)'
u = @. tanh(x^4 - 16*(2*x^2-x'^2))
D = Tridiagonal(ones(m-1),-2ones(m),ones(m-1))/Δx^2
D[1,2] *= 2; D[end,end-1] *= 2
A = I + 0.5Δt*D
B = I - 0.5Δt*D
f = (u,Δt) -> @. u/sqrt(u^2 - (u^2-1)*exp(-50*Δt))
u = f(u,Δt/2)
anim = Animation()
for i = 1:n
(i%10)==1 && (plot(Gray.(u),border=:none); frame(anim))
u = (B\(A*(B\(A*u))'))'
(i<n) && (u = f(u,Δt))
end
u = f(u,Δt/2); Gray.(u)
gif(anim, "allencahn.gif", fps = 30)
14.7. Burgers' equation.
U = (x,t) -> @. t<2 ?
(0<x<t)*(x/t) + (t<x<1+t/2) : (x/t)*(0<x<sqrt(2t))
m = 80; x = LinRange(-1,3,m); Δx = x[2]-x[1]; j = 1:m-1
n = 100; Lₜ = 4; Δt = Lₜ/n
f = u -> u.^2/2; df = u -> u
u = (x.>=0).*(x.<=1)
anim = @animate for i = 1:n
α = 0.5*max.(abs.(df(u[j])),abs.(df(u[j.+1])))
F = (f(u[j])+f(u[j.+1]))/2 - α.*(u[j.+1]-u[j])
global u -= Δt/Δx*[0;diff(F);0]
plot(x,u, fill = (0, 0.3, :blue))
plot!(x,U(x,i*Δt), fill = (0, 0.3, :red))
plot!(legend=:none, ylim=[0,1])
end
gif(anim, "burgers.gif", fps = 15)
14.8. Dam break problem.
δ = u -> diff(u,dims=1)
ϕ = t -> (abs(t)+t)./(1+abs(t))
fixnan(u) = isnan(u)||isinf(u) ? 0 : u
θ = δu -> fixnan.(δu[1:end-1,:]./δu[2:end,:])
∂ₓ(u) = (δu=δ(u);[[0 0];δu[2:end,:].*ϕ.(θ(δu));[0 0]])
F = u -> [u[:,1].*u[:,2] u[:,1].*u[:,2].^2+0.5u[:,1].^2]
m = 100; x = LinRange(-.5,.5,m); Δx = x[2]-x[1]
T = 0.4; n = ceil(T/(Δx/2)); Δt = (T/n)/2;
U = [0.8*(x.<0).+0.2 zero(x)]
Uⱼ = view(U,1:m-1,:); Uⱼ₊₁ = view(U,2:m,:)
anim = @animate for i = 1:n
U° = U-0.5*Δt/Δx*∂ₓ(F(U))
Uⱼ₊₁ .= (Uⱼ+Uⱼ₊₁)/2 - δ(∂ₓ(U))/8 - Δt/Δx*δ(F(U°))
U° = U-0.5*Δt/Δx*∂ₓ(F(U))
Uⱼ .= (Uⱼ+Uⱼ₊₁)/2 - δ(∂ₓ(U))/8 - Δt/Δx*δ(F(U°))
plot(x,U[:,1], fill = (0, 0.3, :blue))
plot!(legend=:none, ylim=[0,1])
end
gif(anim, "dambreak.gif", fps = 15)
15.1. Finite element method.
m=10; x=LinRange(0,1,m); h=x[2]-x[1]
α = fill(2/h-2h/3,m); α[[1,m]] /= 2; β = fill(-1/h-h/6,m-1)
A = SymTridiagonal(α,β)
b = [-2h^3/3; -4h^3/3 .- 8h*x[2:m-1].^2; -4h+8h^2/3-2h^3/3+1]
u = A\b
s = -16 .+ 8x.^2 .+ 15csc(1).*cos.(x)
plot(x,s,marker=:o,alpha=0.5); plot!(x,u,marker=:o,alpha=0.5)
15.2. Finite element method.
m = 8; x = LinRange(0,1,m+2); h = x[2]-x[1]
σ = (a,b,c) -> Tridiagonal(fill(a,m-1),fill(b,m),fill(c,m-1))/h^3
M = [σ(-12,24,-12) σ(-6,0,6);σ(6,0,-6) σ(2,8,2)];
b = [384h*ones(m);zeros(m)]
u = M\b
s = 16*(x.^4 - 2x.^3 + x.^2)
plot(x,[s [0;u[1:m];0]],marker=:o,alpha=0.5)
s = 16*(x.^4 - 2x.^3 + x.^2)
plot(x,[s [0;u[1:m];0]],marker=:o,alpha=0.5)
16.2. Burgers' equation. Fourier spectral methods perform poorly on problems that develop discontinuities such as Burgers' equation. Gibbs oscillations develop around the discontinuity, and these oscillations will spread and grow because Burgers' equation is dispersive. Ultimately, the oscillations overwhelm the solution.
using FFTW, DifferentialEquations
m = 128; x = (1:m)/m*(2π) .- π
ξ = im*[0:(m÷2); (-m÷2+1):-1]
f = (u,p,t) -> -real(ifft(ξ.*fft(0.5u.^2)))
u = solve(ODEProblem(f,exp.(-x.^2),(0.0,2.0)),DP5());
@gif for t=0:.02:2
plot(x,u(t),ylim=(-0.2,1.2))
end
16.3. KdV equation. We'll solve the KdV equation using integrating factors. We first set the initial conditions and parameters. Then, we define the integrating factor G
and the right-hand side f
of the differential equation. Finally, we animate the solution. Notice the two soliton solution.
using FFTW, DifferentialEquations
ϕ = (x,x₀,c) -> 0.5c*sech(sqrt(c)/2*(x-x₀))^2
L = 30; T = 2.0; m = 256
x = (1:m)/m*L .- L/2
u₀ = ϕ.(x,-4,4) + ϕ.(x,-9,9)
iξ = im*[0:(m÷2);(-m÷2+1):-1]*2π/L
F = plan_fft(u₀)
G = t -> exp.(-iξ.^3*t)
f = (w,_,t) -> -G(t) .\ (3iξ .* (F * (F \ (G(t).*w) ).^2) )
w = solve(ODEProblem(f,F*u₀,(0,T)),DP5())
u = t -> real(F\(w(t).*G(t)))
@gif for t=0:.01:2
plot(x,u(t),ylim=(0,5),legend=:none)
end
16.4. Swift–Hohenberg equation. We'll use Strang splitting to solve the Swift–Hohenberg equation.
using FFTW, Images
ϵ = 1; m = 256; ℓ = 100; n = 2000; Δt=100/n
U = (rand(m,m).>.5) .- 0.5
ξ = [0:(m÷2);(-m÷2+1):-1]*2π/ℓ
D²= -ξ.^2 .- (ξ.^2)'
E = exp.(-(D².+1).^2*Δt)
F = plan_fft(U)
f = U -> U./sqrt.(U.^2/ϵ + exp(-Δt*ϵ)*(1 .- U.^2/ϵ))
for i=1:600
U = f(F\(E.*(F*f(U))))
end
Gray.(clamp01.((real(U).+1)/2))
16.5. Cahn–Hilliard equation. We'll solve the two-dimensional Cahn–Hilliard equation using a fourth-order ETDRK scheme.
using FFTW, DifferentialEquations, SciMLOperators
ℓ = 8; T = 8.0; m = 256; ε² = 0.01; α = 1
ξ = [0:(m÷2);(-m÷2+1):-1]*2π/ℓ
Δ = -(ξ.^2 .+ ξ'.^2); Δ² = Δ.^2
u₀ = randn(m,m)
F = plan_fft(u₀)
L = DiagonalOperator((-ε²*Δ²+α*Δ)[:])
N = (u,_,_) -> (v = F\reshape(u,m,m); Δ.*(F*(v.^3-(1+α)*v)))[:]
problem = SplitODEProblem(L,N,(F*u₀)[:],(0,T))
solution = solve(problem,ETDRK4(),dt = 0.1)
u = t -> real(F\reshape(solution(t),m,m))
using Interact, ColorSchemes, Images
@manipulate for t in slider(0:0.1:8; value=0, label="time")
get(colorschemes[:binary], u(t), :extrema)
end
16.6. Kuramoto–Sivashinsky equation. We'll use a fourth-order ETDRK method to solve the two-dimensional KSE.
using FFTW, DifferentialEquations, SciMLOperators
ℓ = 50; T = 200.0; n = 128; u₀ = randn(n,n)
x = LinRange(-ℓ/2,ℓ/2,n+1)[2:end]
ξ = [0:(n÷2);(-n÷2+1):-1]*2π/ℓ
Δ = -(ξ.^2 .+ ξ'.^2); Δ² = Δ.^2
F = plan_fft(u₀)
L = DiagonalOperator((-Δ-Δ²)[:])
N = (u,_,_) ->
( v = reshape(u,n,n);
v = -0.5*abs.((F\(im*ξ.*v)).^2 + (F\(im*ξ'.*v)).^2);
w = (F*v)[:]; w[1] = 0.0; return w )
problem = SplitODEProblem(L,N,(F*u₀)[:],(0,T))
solution = solve(problem, ETDRK4(), dt=0.05, saveat=0.5);
u = t -> real(F\reshape(solution(t),n,n))
using Interact, ColorSchemes, Images
@manipulate for t in slider(0:0.5:T; value=0, label="time")
imresize(get(colorschemes[:magma], -u(t), :extrema),(256,256))
end