{"cells":[{"cell_type": "markdown", "metadata": {}, "source": ["

Numerical Methods in Julia

\n\nThis Jupyter notebook serves as supplementary material to the Julia code from the book [Numerical Methods for Scientific Computing](https://www.equalsharepress.com/media/NMFSC.pdf). 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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using LinearAlgebra, Plots\n","bucket = \"https://raw.githubusercontent.com/nmfsc/data/master/\";"]},{"cell_type": "markdown", "metadata": {}, "source": ["The notebook was written and tested using Julia 1.9.2. Several specialized packages are used throughout this notebook."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["packages = [\"LinearAlgebra\",\"Plots\",\"Images\",\"DelimitedFiles\",\"SparseArrays\",\n"," \"Graphs\",\"GraphPlot\",\"FFTW\",\"Primes\",\"BenchmarkTools\",\"DifferentialEquations\",\n"," \"Optim\",\"Dierckx\",\"Interpolations\",\"SpecialFunctions\",\"Roots\",\"OffsetArrays\",\n"," \"Wavelets\",\"Interact\",\"LsqFit\",\"GLM\",\"DataFrames\",\"Flux\",\"FastGaussQuadrature\",\n"," \"Sundials\",\"SciMLOperators\",\"CSV\",\"ColorSchemes\",\"JuMP\",\"GLPK\",\"MLDatasets\",\n"," \"Arpack\",\"ProgressMeter\",\"IterativeSolvers\",\"Polynomials\",\"ParameterizedFunctions\",\n"," \"ImageShow\",\"NLsolve\",\"DiffEqDevTools\",\"LinRegOutliers\",\"ModelingToolkit\",\"NonlinearSolve\"];\n","#import Pkg; Pkg.add(packages)"]},{"cell_type": "markdown", "metadata": {}, "source": ["

Notebook Contents

\n\n [Part 1: Numerical Linear Algebra](#label0)
\n  [Chapter 1: A Review of Linear Algebra](#label1)
\n  [Chapter 2: Direct Methods for Linear Systems](#label2)
\n  [Chapter 3: Inconsistent Systems](#label3)
\n  [Chapter 4: Computing Eigenvalues](#label4)
\n  [Chapter 6: Fast Fourier Transform](#label5)
\n [Part 2: Numerical Methods for Analysis](#label6)
\n  [Chapter 7: Preliminaries](#label7)
\n  [Chapter 8: Solutions to Nonlinear Equations](#label8)
\n  [Chapter 9: Interpolation](#label9)
\n  [Chapter 10: Approximating Functions](#label10)
\n  [Chapter 11: Differentiation and Integration](#label11)
\n [Part 3: Numerical Differential Equations](#label12)
\n  [Chapter 12: Ordinary Differential Equations](#label13)
\n  [Chapter 13: Parabolic Equations](#label14)
\n  [Chapter 16: Fourier Spectral Methods](#label15)
\n [Part 4: Solutions](#label16)
\n  [Numerical Linear Algebra](#label17)
\n  [Numerical Analysis](#label18)
\n  [Numerical Differential Equations](#label19)
\n"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n# Part 1: Numerical Linear Algebra\n\n## Chapter 1: A Review of Linear Algebra\n**The Hilbert matrix.** The Hilbert matrix $\\mathbf{H}$ is ill-conditioned even for relatively small dimensions. Taking $\\mathbf{H}^{-1}\\mathbf{H}$ should give us the identity matrix. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["hilbert(n) = [1/(i+j-1) for i=1:n, j=1:n]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Images\n","[Gray.(1 .- abs.(hilbert(n)\\hilbert(n))) for n ∈ (10,15,20,25,50)]"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 2: Direct Methods for Linear Systems\n**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. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function gaussian_elimination(A,b)\n"," n = size(A,1)\n"," for j in 1:n\n"," A[j+1:n,j] /= A[j,j]\n"," A[j+1:n,j+1:n] -= A[j+1:n,j:j].*A[j:j,j+1:n]\n"," end\n"," for i in 2:n\n"," b[i:i] -= A[i:i,1:i-1]*b[1:i-1]\n"," end\n"," for i in n:-1:1\n"," b[i:i] = ( b[i] .- A[i:i,i+1:n]*b[i+1:n] )/A[i,i]\n"," end\n"," return b\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = rand(8,8); b = rand(8,1)\n","[A\\b gaussian_elimination(A,b)]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϵ = 1e-15; A = [1 1 2;2 2+ϵ 0; 4 14 4]; b = [-5;10;0.0]\n","A\\b, gaussian_elimination(A,b)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function row_reduce!(tableau)\n"," (i,j) = get_pivot(tableau)\n"," G = tableau[i:i,:]/tableau[i,j]\n"," tableau[:,:] -= tableau[:,j:j]*G\n"," tableau[i,:] = G\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function get_pivot(tableau)\n"," j = argmax(tableau[end,1:end-1])\n"," a, b = tableau[1:end-1,j], tableau[1:end-1,end]\n"," k = findall(a.>0)\n"," i = k[argmin(b[k]./a[k])]\n"," return(i,j) \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function simplex(c,A,b)\n"," (m,n) = size(A)\n"," tableau = [[A I b] ; [c' zeros(1,m) 0]]\n"," while (any(tableau[end,1:n].>0))\n"," row_reduce!(tableau)\n"," end\n"," p = findall(tableau[end,1:n].==0)\n"," x = zeros(n,1)\n"," [x[i]=tableau[:,i]⋅tableau[:,end] for i∈p]\n"," z = -tableau[end,end]\n"," y = -tableau[end,n.+(1:m)]\n"," return((z=z, x=x, y=y))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = [2 0 1;4 1 2;1 1 0]; c = [2;1;1]; b = [3;2;1]\n","solution = simplex(c,A,b)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Graph drawing.** The following code draws the dolphin networks of the Doubtful Sound. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DelimitedFiles, SparseArrays\n","function get_adjacency_matrix(filename)\n"," ij = readdlm(download(bucket*filename*\".csv\"),',',Int)\n"," sparse(ij[:,1],ij[:,2],one.(ij[:,1]))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Graphs, GraphPlot\n","g = get_adjacency_matrix(\"dolphins\") |> Graph\n","gplot(g,layout=spectral_layout)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["gplot(g,layout=spring_layout)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["gplot(g,layout=circular_layout)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Revised simplex method.** "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using SparseArrays\n","function revised_simplex(c,A,b)\n"," (m,n) = size(A)\n"," ξ = i -> (z=spzeros(m);z[i]=1;z)\n"," N = Vector(1:n); B = Vector(n .+ (1:m))\n"," A = [A sparse(I, m, m)]\n"," ABinv = sparse(I, m, m)\n"," c = [c;spzeros(m,1)]\n"," while(true)\n"," j = findfirst(x->x>0,(c[N]'.-(c[B]'*ABinv)*A[:,N])[:])\n"," if isnothing(j); break; end\n"," q = ABinv*A[:,N[j]]\n"," k = findall(q.>0)\n"," i = k[argmin(ABinv[k,:]*b./q[k])]\n"," B[i], N[j] = N[j], B[i]\n"," ABinv -= ((q - ξ(i))/q[i])*ABinv[i:i,:]\n"," end\n"," i = findall(B.≤n)\n"," x = zeros(n,1)\n"," x[B[i]] = ABinv[i,:]*b\n"," return((z=c[1:n]'*x, x=x, y=c[B]'*ABinv))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = [2 0 1;4 1 2;1 1 0]; c = [2;1;1]; b = [3;2;1]\n","solution = simplex(c,A,b)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 3: Inconsistent Systems\n**Zipf's law.** Let's use ordinary least squares to find Zipf's law coefficients for the canon of Sherlock Holmes."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DelimitedFiles\n","T = readdlm(download(bucket*\"sherlock.csv\"), '\\t')[:,2]\n","n = length(T)\n","A = [ones(n,1) log.(1:n)]\n","B = log.(T)\n","c = A\\B\n","print(\"ordinary least squares:\\n\"*string(c)*\"\\n\")"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["zipf(x,c) = exp(c[1])*x.^c[2]\n","scatter(T,xaxis=:log,yaxis=:log,ma=0.4,msw=0)\n","plot!([1,n],zipf([1,n],c),w=3)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Constrained least squares.** The constrained least squares problem of solving $\\mathbf{Ax} = \\mathbf{b}$ with the constraint condition $\\mathbf{Cx}=\\mathbf{d}$:"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function constrained_lstsq(A,b,C,d)\n"," x = [A'*A C'; C zeros(size(C,1),size(C,1))]\\[A'*b;d]\n"," x[1:size(A,2)]\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Total least squares.** The function `tls` solves the total least squares problem."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function tls(A,B)\n"," n = size(A,2)\n"," V = svd([A B]).V\n"," -V[1:n,n+1:end]/V[n+1:end,n+1:end]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = [2 4; 1 -1; 3 1; 4 -8]; b = [1; 1; 4; 1];\n","xₒₗₛ = A\\b; xₜₗₛ = tls(A,b)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["print(\"ordinary least squares:\"*string(xₒₗₛ))\n","print(\"\\ntotal least squares:\"*string(xₜₗₛ))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Images\n","A = load(download(bucket*\"laura.png\")) .|> Gray\n","U, σ, V = svd(A);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["k = 20\n","Aₖ = U[:,1:k] * Diagonal(σ[1:k]) * V[:,1:k]' .|> Gray\n","[A Aₖ]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["norm(A-Aₖ) ≈ norm(σ[k+1:end])"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϵ² = 1 .- cumsum(σ)/sum(σ); scatter(ϵ²,xaxis=:log)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Nonnegative matrix factorization.** The following code is a naive implementation of nonnegative matrix factorization using multiplicative updates (without a stopping criterion):"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function nmf(X,p=6)\n"," W = rand(Float64, (size(X,1), p))\n"," H = rand(Float64, (p,size(X,2)))\n"," for i in 1:50\n"," W = W.*(X*H')./(W*(H*H') .+ (W.≈0))\n"," H = H.*(W'*X)./((W'*W)*H .+ (H.≈0))\n"," end\n"," (W,H)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Images\n","A = Gray.(load(download(bucket*\"laura.png\")))\n","W,H = nmf(Float64.(A),20)\n","[A Gray.(W*H)]"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 4: Computing Eigenvalues\n**Eigenvalue condition number.** The function `condeig` computes the eigenvalue condition number. Let's use it on a small random matrix."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function condeig(A)\n"," Sᵣ = eigvecs(A)\n"," Sₗ = inv(Sᵣ')\n"," Sₗ ./= sqrt.(sum(abs.(Sₗ.^2), dims=1))\n"," 1 ./ abs.(sum(Sᵣ.*Sₗ, dims=1))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["condeig(rand(4,4))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**PageRank.** The following minimal code computes the PageRank of the very small graph by using the power method over 9 iterations
\"internet"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["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]\n","v = all(H.==0,dims=1)\n","H = H ./ (sum(H,dims=1)+v) \n","n = size(H,1) \n","d = 0.85\n","x = ones(n,1)/n\n","for i in 1:9\n"," x = d*(H*x) .+ d/n*(v*x) .+ (1-d)/n\n","end \n","x"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, Images, ColorSchemes\n","perfectshuffle= A -> [A[1:2:end,:];A[2:2:end,:]]\n","F = real.(fft(I(128),1)); PF = F\n","for i = 1:6\n"," F = perfectshuffle(F); PF = [PF F]\n","end\n","get(colorschemes[:gray1], PF, :extrema)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 6: Fast Fourier Transform\n**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`."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function fftx2(c)\n"," n = length(c)\n"," ω = exp(-2im*π/n)\n"," if mod(n,2) == 0 \n"," k = collect(0:n/2-1)\n"," u = fftx2(c[1:2:n-1])\n"," v = (ω.^k).*fftx2(c[2:2:n]) \n"," return([u+v; u-v])\n"," else \n"," k = collect(0:n-1)\n"," F = ω.^(k*k')\n"," return(F*c)\n"," end \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ifftx2(y) = conj(fftx2(conj(y)))/length(y);"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Fast Toeplitz multiplication.** The function `fasttoeplitz` computes the Toeplitz multiplication by padding out a Toeplitz matrix with zeros to make it circulant."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function fasttoeplitz(c,r,x)\n"," n = length(x)\n"," Δ = nextpow(2,n) - n\n"," x₁ = [c; zeros(Δ); r[end:-1:2]]\n"," x₂ = [x; zeros(Δ+n-1)]\n"," ifftx2(fftx2(x₁).*fftx2(x₂))[1:n]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, Primes, Plots\n","N = 10000\n","smooth(n,N) = (1:N)[all.(x->x<=n,factor.(Set,1:N))]\n","t₁ = [(x = randn(n); (n,@elapsed(fft(x)))) for n∈primes(N)]\n","t₂ = [(x = randn(n); (n,@elapsed(fft(x)))) for n∈smooth(5,N)]\n","plot(t₁,label=\"prime\"); plot!(t₂,label=\"5-smooth\")\n","plot!(ylims=(0,t₁[end][2]*3))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function bluestein_fft(x)\n"," n = length(x)\n"," ω = exp.((1im*π/n)*(0:n-1).^2)\n"," ω.*fasttoeplitz(conj(ω),conj(ω),ω.*x)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Primes\n","function primitiveroot(n)\n"," ϕ = n - 1\n"," p = factor(Set,ϕ)\n"," for r = 2:n-1\n"," all( [powermod(r, ϕ÷pᵢ, n) for pᵢ∈p] .!= 1 ) && return(r)\n"," end\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function rader_fft(x)\n"," n = length(x)\n"," r = primitiveroot(n)\n"," P₊ = powermod.(r, 0:n-2, n)\n"," P₋ = circshift(reverse(P₊),1)\n"," ω = exp.((2im*π/n)*P₋)\n"," c = x[1] .+ ifft(fft(ω).*fft(x[2:end][P₊]))\n"," [sum(x); c[reverse(invperm(P₋))]] \n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW\n","dst(x) = FFTW.r2r(x,FFTW.RODFT00)\n","idst(x) = dst(x)/(2^ndims(x)*prod(size(x).+1))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 50; x = (1:n)/(n+1); Δx = 1/(n+1) \n","v = 2 .- 2cos.(x*π) \n","λ = [v[i]+v[j]+v[k] for i∈1:n, j∈1:n, k∈1:n]./Δx^2\n","f = [(x-x^2)*(y-y^2) + (x-x^2)*(z-z^2)+(y-y^2)*(z-z^2) \n"," for x∈x,y∈x,z∈x]\n","u = idst(dst(f)./λ);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["norm(u - [(x-x^2)*(y-y^2)*(z-z^2)/2 for x∈x,y∈x,z∈x])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Image compression.** Let's consider the following code that compresses a grayscale image `A` to a factor of `d` from its original size. "]},{"cell_type": "markdown", "metadata": {}, "source": ["**DCT image compression.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Images, FFTW, SparseArrays\n","function dctcompress(A,d)\n"," B = dct(Float64.(A))\n"," idx = d*prod(size(A)) |> floor |> Int\n"," tol = sort!(abs.(B[:]),rev=true)[idx]\n"," B₀ = droptol!(sparse(B),tol)\n"," idct(B₀) |> clamp01! .|> Gray\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = Gray.(load(download(bucket*\"laura_square.png\")))\n","[A dctcompress(A,0.05)]"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n# Part 2: Numerical Methods for Analysis\n\n## Chapter 7: Preliminaries\nLet's start with a function that returns a double-precision floating-point representation as a string of bits."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["bitstring(Float64(π))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Fast inverse square root.** The following function implements the circa 1999 Q_rsqrt algorithm to approximate the reciprocal square root of a number."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function Q_rsqrt(x::Float32)\n"," i = reinterpret(Int32,Float32(x)) \n"," i = Int32(0x5f3759df) - (i>>1) \n"," y = reinterpret(Float32,i) \n"," y = y * (1.5f0 - (0.5f0 * x * y * y)) \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Q_rsqrt(Float32(0.01))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using BenchmarkTools\n","@btime Q_rsqrt(Float32(x)) setup=(x=rand()) seconds=3\n","@btime 1/sqrt(x) setup=(x=rand()) seconds=3"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Rump's catastrophic cancellation.** The answer should be `-0.827396...` What does Julia come up with?"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["a = 77617; b = 33096\n","333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b)"]},{"cell_type": "markdown", "metadata": {}, "source": ["`NaN` can be used to lift \"pen off paper\" when plotting a series of connected points."]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 8: Solutions to Nonlinear Equations\nWe start with simple implementation of the bisection method."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function bisection(f,a,b,tolerance)\n"," while abs(b-a) > tolerance\n"," c = (a+b)/2\n"," sign(f(c)) == sign(f(a)) ? a = c : b = c\n"," end\n"," (a+b)/2\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["bisection(x->sin(x),2,4,1e-14)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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\\}$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function escape(n, c=0, z=0)\n"," for k in 1:n\n"," z = z^2 + c\n"," abs2(z) > 4 && return k\n"," end\n"," return n\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function mandelbrot(bb, xn=800, n=200, z=0)\n"," yn = (xn*(bb[4]-bb[2]))÷(bb[3]-bb[1]) |> Int\n"," c = LinRange(bb[1],bb[3],xn)' .+ im*LinRange(bb[4],bb[2],yn)\n"," [escape(n,c,z) for c∈c]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Images\n","M = mandelbrot([-0.172, 1.0228, -0.1494, 1.0443])\n","save(\"mandelbrot.png\",1 .- M./maximum(M))\n","load(\"mandelbrot.png\")"]},{"cell_type": "markdown", "metadata": {}, "source": ["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`."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function julia(bb, xn=800, n=200, c=-0.123+0.745im)\n"," yn = (xn*(bb[4]-bb[2]))÷(bb[3]-bb[1]) |> Int\n"," z = LinRange(bb[1],bb[3],xn)' .+ im*LinRange(bb[4],bb[2],yn)\n"," [escape(n,c,z) for z∈z]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["J = julia([-2, -2, 2, 2])\n","save(\"julia.png\",1 .- J./maximum(J))\n","load(\"julia.png\")"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Homotopy continuation.** The following snippet of code finds a root of $$x^3-3xy^2-1 =0$$\n$$y^3-3x^2y = 0$$ with an initial guess $(x,y) = (1,1)$ using homotopy continuation: "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations\n","f = z -> ((x,y)=tuple(z...); [x^3-3x*y^2-1; y^3-3x^2*y] )\n","df = (z,p,_) -> ((x,y)=tuple(z...);\n"," -[3x^2-3y^2 -6x*y; -6x*y 3y^2-3x^2]\\p)\n","z₀ = [1,1]\n","sol = solve(ODEProblem(df,z₀,(0,1),f(z₀)))\n","sol.u[end]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function plot_trace(x,f)\n"," contour(-2:0.02:1, -0.25:0.02:3.25,(x,y) -> cbrt(f([x,y])), color=:red, colorbar=:none)\n"," plot!(x[1,:],x[2,:],color=:black, linewidth=2, marker=:o, legend=:none)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["∇f = x-> [-2(1-x[1])-400x[1]*(x[2]-x[1]^2); 200(x[2]-x[1]^2)]\n","x = [-1.8,3.0]; p = [0,0]; α = 0.001; β = 0.9\n","s = x\n","for i = 1:100\n"," p = -∇f(x) + β*p\n"," x += α*p\n"," append!(s,x)\n","end\n","f = x -> (1-x[1])^2 + 100(x[2] - x[1]^2)^2\n","plot_trace(reshape(s, 2, :),f)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Optim\n","x₀ = [-1.8,3.0]\n","options = Optim.Options(store_trace = true,extended_trace = true)\n","result = optimize(f, x₀, BFGS(),options)\n","s = Optim.x_trace(result); s = vcat(s...)\n","plot_trace(reshape(s, 2, :),f)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["∇f = x -> [-2(1-x[1])-400x[1]*(x[2]-x[1]^2); 200(x[2]-x[1]^2)]\n","x = [-1.8,3.0]; α = 0.001; p = [0,0]; β = 0.9\n","for i = 1:100\n"," p = -∇f(x) + β*p\n"," x += α*p\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Optim\n","f = x -> (1-x[1])^2 + 100(x[2] - x[1]^2)^2\n","x₀ = [-1.8,3.0]\n","result = optimize(f, x₀, GradientDescent())\n","x = Optim.minimizer(result)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 9: Interpolation\n **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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function spline_natural(xᵢ,yᵢ)\n"," h = diff(xᵢ)\n"," γ = 6*diff(diff(yᵢ)./h)\n"," α = h[2:end-1]\n"," β = 2(h[1:end-1]+h[2:end])\n"," [0;SymTridiagonal(β,α)\\γ;0]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function evaluate_spline(xᵢ,yᵢ,m,n)\n"," h = diff(xᵢ)\n"," B = yᵢ[1:end-1] .- m[1:end-1].*h.^2/6\n"," A = diff(yᵢ)./h - h./6 .*diff(m)\n"," x = range(minimum(xᵢ),maximum(xᵢ),length=n+1)\n"," i = sum(xᵢ'.≤x,dims=2)\n"," i[end] = length(xᵢ)-1 \n"," y = @. (m[i]*(xᵢ[i+1]-x)^3 + m[i+1]*(x-xᵢ[i])^3)/(6h[i]) +\n"," A[i]*(x-xᵢ[i]) + B[i]\n"," return(x,y)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["x = LinRange(0,1,8); y = rand(8)\n","m = spline_natural(x,y)\n","X,Y = evaluate_spline(x,y,m,100)\n","scatter(x,y); plot!(X,Y)"]},{"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["g = x-> @. max(1-abs(3-x),0)\n","xᵢ = 0:5; yᵢ = g(xᵢ)\n","x = LinRange(0,5,101);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Dierckx, Plots\n","spline = Spline1D(xᵢ,yᵢ)\n","plot(x,spline(x)); plot!(x,g(x)); scatter!(xᵢ,yᵢ)"]},{"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Interpolations\n","method = BSpline(Cubic(Natural(OnGrid())))\n","spline = scale(interpolate(yᵢ, method), xᵢ)\n","plot(x,spline(x)); plot!(x,g(x)); scatter!(xᵢ,yᵢ)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["bernstein(n,t) = @. binomial(n,0:n)'*t^(0:n)'*(1-t)^(n:-1:0)'"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 3\n","t = LinRange(0,1,100)\n","p = rand(n+1,2)\n","z = bernstein(n,t)*p\n","plot(p[:,1],p[:,2],marker=:o,opacity=0.3); \n","plot!(z[:,1],z[:,2],legend=:none)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 20; f = t -> sqrt(t)\n","y = bernstein(n,t)*f.(LinRange(0,1,n+1))\n","plot(t,y); plot!(t,f.(t))"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 10: Approximating Functions\n **Legendre polynomials.** We can evaluate a Legendre polynomial of order $n$ using Bonnet's recursion formula."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function legendre(x,n)\n"," n==0 && return(one.(x))\n"," n==1 && return(x)\n"," x.*legendre(x,n-1) .- (n-1)^2/(4(n-1)^2-1)*legendre(x,n-2)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["x = LinRange(-1,1,100)\n","plot()\n","for n in 0:4\n"," plot!(x,legendre(x,n)) \n","end\n","current()"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Chebyshev polynomials.** We'll construct a Chebyshev differentiation matrix and use the matrix to solve a few simple problems."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function chebdiff(n)\n"," x = -cos.(LinRange(0,π,n))\n"," c = [2;ones(n-2);2].*(-1).^(0:n-1)\n"," D = c./c'./(x .- x' + I)\n"," D - Diagonal([sum(D,dims=2)...]), x \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 15\n","D,x = chebdiff(n)\n","u = exp.(-4x.^2)\n","plot(x,D*u,marker=:o)\n","t = LinRange(-1,1,200)\n","plot!(t,-8t.*exp.(-4t.^2))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["B = zeros(1,n); B[1] = 1\n","plot(x,[B;D]\\[2;u],marker=:o)"]},{"cell_type": "markdown", "metadata": {}, "source": ["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$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 15; k² = 256\n","D,x = chebdiff(n)\n","L = (D^2 - k²*Diagonal(x))\n","L[[1,n],:] .= 0; L[1,1] = L[n,n] = 1 \n","y = L\\[2;zeros(n-2);1]\n","plot(x,y,marker=:o)"]},{"cell_type": "markdown", "metadata": {}, "source": ["The analytical solution can be expressed as the sum of Airy functions."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using SpecialFunctions\n","a = [airyai(-∛k²) airybi(-∛k²);airyai(∛k²) airybi(∛k²)]\\[2;1]\n","airy(x) = a[1]*airyai(∛k²*x) + a[2]*airybi(∛k²*x)\n","scatter(x,y)\n","plot!(t,airy.(t))"]},{"cell_type": "markdown", "metadata": {}, "source": ["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?"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["λ² = 16^2; N = 6:2:60; ϵ = []\n","for n ∈ N \n"," D,x = chebdiff(n)\n"," L = (D^2 - k²*Diagonal(x))\n"," L[[1,n],:] .= 0; L[1,1] = L[n,n] = 1 \n"," y = L\\[2;zeros(n-2);1]\n"," append!(ϵ,norm(y - airy.(x),Inf))\n","end\n","plot(N,ϵ,yaxis=:log,marker=:o,legend=:none)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using OffsetArrays\n","function scaling(c,ϕₖ,n)\n"," m = length(c)-1; ℓ = 2^n\n"," ϕ = OffsetVector(zeros(3*m*ℓ),-m*ℓ) \n"," k = (0:m)*ℓ\n"," ϕ[k] = ϕₖ\n"," for j = 1:n\n"," for i = 1:m*2^(j-1)\n"," x = (2i-1)*2^(n-j) \n"," ϕ[x] = c ⋅ ϕ[2x .- k]\n"," end\n"," end\n"," ((0:m*ℓ-1)/ℓ, ϕ[0:m*ℓ-1])\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["c = [1+√3, 3+√3, 3-√3, 1-√3]/4\n","z = [0, 1+√3, 1-√3, 0]/2\n","(x,ϕ) = scaling(c, z, 8)\n","plot(x,ϕ)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ψ = zero(ϕ); n = length(c)-1; ℓ = length(ψ)÷2n\n","for k∈0:n\n"," ψ[(k*ℓ+1):(k+n)*ℓ] += (-1)^k*c[n-k+1]*ϕ[1:2:end]\n","end\n","plot(x,ψ)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**DWT image compression.** The following code explores using a discrete wavelet transform along with filtering as means of image compression."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Wavelets, Images\n","img = float.(Gray.(load(download(bucket*\"laura_square.png\"))))\n","B = dwt(Float32.(img), wavelet(WT.haar))\n","[img Gray.(1 .- clamp01.(sqrt.(abs.(B))))]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function dwtfilter(channels,wt,k) \n"," map(1:3) do i\n"," A = dwt(channels[i,:,:], wavelet(wt))\n"," threshold!(A, HardTH(), k) \n"," clamp01!(idwt(A, wavelet(wt)))\n"," end\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Interact, Wavelets, Images\n","img = float.(load(download(bucket*\"laura_square.png\")))\n","channels = channelview(img)\n","func = Dict(\"Haar\"=>WT.haar, \"D₄\"=>WT.db4, \"Coiflet\"=>WT.coif4)\n","@manipulate for wt in togglebuttons(func; label=\"Transform\"),\n"," k in slider(0:0.05:1.0; value = 0, label = \"Threshold\")\n"," colorview(RGB,dwtfilter(channels,wt,k)...)\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function jacobian(f,x,c)\n"," J = zeros(length(x),length(c))\n"," for k in (n = 1:length(c))\n"," J[:,k] .= imag(f(x,c+1e-8im*(k .== n)))/1e-8\n"," end\n"," return J\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function gauss_newton(f,x,y,c)\n"," r = y - f(x,c) \n"," for j = 1:100\n"," G = jacobian(f,x,c)\n"," M = G'*G\n"," c += (M+Diagonal(M))\\(G'*r)\n"," norm(r-(r=y-f(x,c))) < 1e-12 && return(c)\n"," end\n"," print(\"Gauss-Newton did not converge.\")\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["f = (x,c) -> @. c[1]*exp(-c[2]*(x-c[3])^2) +\n"," c[4]*exp(-c[5]*(x-c[6])^2)\n","x = 8*rand(100)\n","y = f(x,[1, 3, 3, 2, 3, 6]) + 0.1*randn(100);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["c₀ = [2, 0.3, 2, 1, 0.3, 7]\n","c = gauss_newton(f,x,y,c₀)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["if !isnothing(c)\n"," X = LinRange(0,8,200)\n"," scatter(x,y,opacity=0.5)\n"," plot!(X,f(X,c))\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["In practice, we might use the LsqFit.jl package which solves nonlinear least squares problems."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using LsqFit\n","cf = curve_fit(f,x,y,c₀)\n","c, r = cf.param, cf.resid\n","scatter(x,y,opacity=0.5)\n","plot!(X,f(X,c))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["σ = x -> @. 1/(1+exp(-x))\n","x = rand(30); y = ( rand(30) .< σ(10x.-5) );"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["X, w = [one.(x) x], zeros(2,1)\n","for i=1:10\n"," S = σ(X*w).*(1 .- σ(X*w))\n"," w += (X'*(S.*X))\\(X'*(y - σ(X*w)))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["scatter(x,y)\n","t = LinRange(minimum(x),maximum(x),50)\n","plot!(t,σ([one.(t) t]*w))"]},{"cell_type": "markdown", "metadata": {}, "source": ["In practice, we might use Julia's GLM library to solve the logistic regression problem."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using GLM, DataFrames\n","data = DataFrame(X=x, Y=y)\n","logit = glm(@formula(Y ~ X), data, Bernoulli(), LogitLink())"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Neural Networks.** Let's use a neural network to find a function that approximates semi-circular data."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["N = 100; θ = LinRange(0,π,N)'\n","x = cos.(θ); x̃ = [one.(x);x]\n","y = sin.(θ) + 0.05*randn(1,N)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 20; W₁ = rand(n,2); W₂ = randn(1,n)\n","ϕ = x -> max.(x,0)\n","dϕ = x -> (x.>0)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["α = 0.1\n","for epoch = 1:1000\n"," ŷ = W₂ * ϕ(W₁*x̃)\n"," ∂L∂y = 2(ŷ-y)/N\n"," ∂L∂W₁ = dϕ(W₁*x̃) .* (W₂' * ∂L∂y) * x̃'\n"," ∂L∂W₂ = ∂L∂y * ϕ(W₁*x̃)'\n"," W₁ -= 0.1α * ∂L∂W₁ \n"," W₂ -= α * ∂L∂W₂ \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["scatter(x̃[2,:],y',opacity=0.3)\n","x̂ = LinRange(-1.2,1.2,200)'; x̂ = [one.(x̂);x̂]; ŷ = W₂ * ϕ(W₁*x̂)\n","plot!(x̂[2,:],ŷ',width=3)"]},{"cell_type": "markdown", "metadata": {}, "source": ["We can solve the same problem using a sigmoid activation function:"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 20; W₁ = rand(n,2); W₂ = randn(1,n)\n","ϕ = x -> @. 1 / (1 + exp(-x))\n","dϕ = x -> @. ϕ(x)*(1-ϕ(x))\n","α = 0.1\n","for epoch = 1:10000\n"," ŷ = W₂ * ϕ(W₁*x̃)\n"," L = norm(ŷ - y)\n"," ∂L∂y = (ŷ-y)/L\n"," ∂L∂W₁ = dϕ(W₁*x̃) .* (W₂' * ∂L∂y) * x̃'\n"," ∂L∂W₂ = ∂L∂y * ϕ(W₁*x̃)'\n"," W₁ -= 0.1α * ∂L∂W₁ \n"," W₂ -= α * ∂L∂W₂ \n","end\n","scatter(x̃[2,:],y',opacity=0.3)\n","x̂ = LinRange(-1.2,1.2,200)'; x̂ = [one.(x̂);x̂]; ŷ = W₂ * ϕ(W₁*x̂)\n","plot!(x̂[2,:],ŷ',width=3)"]},{"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Flux\n","model = Chain(Dense(1, n, relu), Dense(n, 1))\n","loss(x,y) = Flux.Losses.mse(model(x), y)\n","parameters = Flux.params(model)\n","data = [(x,y)]\n","optimizer = Descent(0.1)\n","for epochs = 1:10000\n"," Flux.train!(loss, parameters, data, optimizer)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["scatter(x',y',ma=0.4,msw=0)\n","plot!(x',model(x)',width=3)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 11: Differentiation and Integration\n Let's compute the coefficients to the third-order approximation to $f'(x)$ using nodes at $x-h$, $x$, $x+h$ and $x+2h$. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["δ = [-1,0,1,2]\n","n = length(δ)\n","V = δ.^(0:n-1)' .// factorial.(0:n-1)'\n","C = inv(V)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Richardson extrapolation.** $D(\\phi(x))$ of a finite difference operator $\\phi(x)$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function richardson(f,x,m,n=m)\n"," n > 0 ? \n"," (4^n*richardson(f,x,m,n-1) - richardson(f,x,m-1,n-1))/(4^n-1) :\n"," ϕ(f,x,2^m) \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϕ = (f,x,n) -> (f(x+1/n) - f(x-1/n))/(2/n)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["richardson(x->sin(x),0,4)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Automatic differentiation.** We can build a minimal working example of forward accumulation automatic differentiation by defining a class and overloading the base operators. We'll verify the code using the function $x + \\sin x$. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["struct Dual\n"," value\n"," deriv\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Dual(x) = Dual(x,1)\n","value(x) = x isa Dual ? x.value : x\n","deriv(x) = x isa Dual ? x.deriv : 0"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Base.:+(u, v) = Dual(value(u)+value(v), deriv(u)+deriv(v))\n","Base.:-(u, v) = Dual(value(u)-value(v), deriv(u)-deriv(v))\n","Base.:*(u, v) = Dual(value(u)*value(v), \n"," value(u)*deriv(v)+value(v)*deriv(u))\n","Base.:sin(u) = Dual(sin(value(u)), cos(value(u))*deriv(u))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["x = Dual(0)\n","y = x + sin(x)"]},{"cell_type": "markdown", "metadata": {}, "source": ["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)$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["x₁ = Dual(2,[1 0])\n","x₂ = Dual(π,[0 1])\n","y₁ = x₁*x₂ + sin(x₂)\n","y₂ = x₁*x₂ - sin(x₂)\n","value(y₁),value(y₁),deriv(y₁),deriv(y₂)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function trapezoidal(f,x,n) \n"," F = f.(LinRange(x[1],x[2],n+1))\n"," (F[1]/2 + sum(F[2:end-1]) + F[end]/2)*(x[2]-x[1])/n\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϕ = (f,x,n) -> trapezoidal(f,x,n)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["richardson(x->sin(x),[0,π/2],4)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = [floor(Int,10^y) for y in LinRange(1, 2, 10)]\n","error = zeros(10,7)\n","f = (x,p) -> x + x.^p.*(2-x).^p\n","for p ∈ 1:7,\n"," S = trapezoidal(x->f(x,p),[0,2],10^6)\n"," for i ∈ 1:length(n)\n"," Sₙ = trapezoidal(x->f(x,p),[0,2],n[i])\n"," error[i,p] = abs(Sₙ - S)/S\n"," end\n","end\n","slope = ([log.(n) one.(n)]\\log.(error))[1,:]\n","info = .*(string.((1:7)'),\": slope=\",string.(round.(slope')))\n","plot(n,error, xaxis=:log, yaxis=:log, labels = info)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, LinearAlgebra\n","function clenshaw_curtis(f,n)\n"," x = cos.(π*(0:n)'/n)\n"," w = zeros(n+1,1); w[1:2:n+1] = 2 ./ (1 .- (0:2:n).^2)\n"," 1/n * dctI(f.(x)) ⋅ w\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function dctI(f) \n"," g = FFTW.r2r!([f...],FFTW.REDFT00)\n"," [g[1]/2; g[2:end-1]; g[end]/2] \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["clenshaw_curtis(x -> cos(8x^2),20)"]},{"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function dctI(f)\n"," n = length(f)\n"," g = real(fft(f[[1:n; n-1:-1:2]]))\n"," [g[1]/2; g[2:n-1]; g[n]/2] \n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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$. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function gauss_legendre(n)\n"," a = zeros(n) \n"," b = (1:n-1).^2 ./ (4*(1:n-1).^2 .- 1)\n"," 𝟙² = 2\n"," λ, v = eigen(SymTridiagonal(a, sqrt.(b))) \n"," (λ, 𝟙²*v[1,:].^2)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 8\n","f = x -> cos(x)*exp(-x^2);\n","nodes, weights = gauss_legendre(n)\n","f.(nodes) ⋅ weights"]},{"cell_type": "markdown", "metadata": {}, "source": ["Alternatively, the Julia library FastGaussQuadrature.jl provides a fast, accurate method of computing the nodes and weights"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FastGaussQuadrature\n","nodes, weights = gausslegendre(n)\n","f.(nodes) ⋅ weights"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n# Part 3: Numerical Differential Equations\n\n## Chapter 12: Ordinary Differential Equations\n Let's plot the boundary of the region of absolute stability for BDF2: $$z = \\frac{\\frac{3}{2} r^2 - 2 r + \\frac{1}{2}}{r^2}$$"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["r = exp.(2im*π*(0:.01:1))\n","plot(@. (1.5r^2 - 2r + 0.5)/r^2); plot!(aspect_ratio=:equal)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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]`."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function multistepcoefficients(m,n)\n"," s = length(m) + length(n) - 1\n"," A = (m.+1).^(0:s-1) .|> Rational\n"," B = (0:s-1).*((n.+1).^[0;0:s-2])\n"," c = -[A[:,2:end] B]\\ones(Int64,s)\n"," [1;c[1:length(m)-1]], c[length(m):end] \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function plotstability(a,b)\n"," λk(r) = (a ⋅ r.^-(0:length(a)-1)) ./ (b ⋅ r.^-(0:length(b)-1))\n"," r = exp.(im*LinRange(0,2π,200))\n"," plot(λk.(r),label=\"\",aspect_ratio=:equal)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m = [0 1]; n = [0 1 2]\n","a, b = zeros(maximum(m)+1), zeros(maximum(n)+1) \n","a[m.+1], b[n.+1] = multistepcoefficients(m,n)\n","plotstability(a,b)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Recipe for solving an ODE.** The general recipe for solving an ODE is to\n1. Load the module\n1. Set up the parameters\n1. Define the problem\n1. Choose the method\n1. Solve the problem\n1. Present the solution\n\nLet'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]$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations, Plots\n","pendulum(u,p,t) = [u[2]; -sin(u[1])]\n","u₀ = [8π/9,0]; tspan = [0,8π]\n","problem = ODEProblem(pendulum, u₀, tspan)\n","method = Trapezoid()\n","solution = solve(problem,method)\n","plot(solution, xaxis=\"t\", label=[\"θ\" \"ω\"])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations, ParameterizedFunctions, DiffEqDevTools, Sundials\n","f = @ode_def LotkaVolterra begin\n"," dx = α*x - β*x*y\n"," dy = -γ*y + δ*x*y\n","end α β γ δ\n","\n","params = [1.5,1.0,3.0,1.0]\n","problem = ODEProblem(f,[1.0;1.0],(0.0,10.0),params)\n","\n","solution = solve(problem,Vern7(),abstol=1e-15,reltol=1e-15)\n","testsol = TestSolution(solution)\n","plot(solution)"]},{"cell_type": "markdown", "metadata": {}, "source": ["Now, we'll choose a set of solvers to benchmark."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["solvers = [ \n"," Dict(:alg=>VCABM())\n"," Dict(:alg=>DP5())\n"," Dict(:alg=>BS3())\n"," Dict(:alg=>CVODE_Adams())\n"," Dict(:alg=>Tsit5())\n"," Dict(:alg=>DP8())\n"," Dict(:alg=>Vern8())\n"," Dict(:alg=>Vern9())\n","];"]},{"cell_type": "markdown", "metadata": {}, "source": ["We use DiffEqDevTools.jl to generate a working precision set."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["abstols = 10 .^[-(LinRange(6,13,40))...]\n","reltols = 1e3*abstols\n","wps = WorkPrecisionSet(problem,abstols,reltols,solvers;\n"," appxsol=testsol,save_everystep=true,numruns=100,maxiters=10^5)"]},{"cell_type": "markdown", "metadata": {}, "source": ["Let's plot the data and fit regression lines. First, we'll remove the outliers."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DataFrames, LinRegOutliers\n","function logfit(w, c = 1)\n"," df = DataFrame(x=log.(w.errors), y=log.(w.times))\n"," reg = createRegressionSetting(@formula(y ~ x), df)\n"," s = lts(reg; crit = c)\n"," β = s[\"betas\"]\n"," delete!(df, s[\"outliers\"])\n"," x = exp.([minimum(df.x) maximum(df.x)]) \n"," y = [exp(β[1])*(x[1])^β[2] exp(β[1])*(x[2])^β[2]] \n"," (x,y,exp.(df.x),exp.(df.y))\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["Now, we plot the compute time versus precision."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["plot(xaxis=:log,yaxis=:log,xflip = true)\n","for i∈1:length(wps)\n"," x, y, xₚ, yₚ = logfit(wps[i],5)\n"," scatter!(xₚ,yₚ,alpha=0.25,labels=:none,markersize=5,color=i)\n"," plot!(x',y',linewidth=3,labels=:none,color=i)\n"," annotate!(x[1], y[1], text.(wps.names[i],9,:left)) \n","end\n","plot!()"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Differential-algebraic equations.** The pendulum problem."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations\n","θ₀ = π/3; ℓ = 1; tspan = (0,30.0)\n","U₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0]\n","function pendulum(dU,U,p,t)\n"," x,y,u,v = U; ℓ,g,m = p\n"," τ = m*(u^2 + v^2 + g*y)/ℓ^2\n"," dU .= (u, v, -τ*x/m, -τ*y/m + g)\n","end\n","problem = ODEProblem(pendulum, U₀, tspan, (ℓ,1,1))\n","solution = solve(problem, Tsit5())\n","plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function residual(r,w,p,t)\n"," x,y,u,v = w; ℓ,g,m = p\n"," r .= (x^2+y^2-ℓ^2, x*u+y*v, 0, 0)\n","end\n","cb = ManifoldProjection(residual)\n","solution = solve(problem, Tsit5(), callback=cb)\n","plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations\n","θ₀ = π/3; ℓ = 1; tspan = (0,30.0)\n","u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, -ℓ*cos(θ₀)]\n","function pendulum(dU,U,p,t)\n"," x,y,u,v,τ = U; ℓ,g,m = p\n"," dU .= (u, v, -τ*x/m, -τ*y/m + g, -ℓ^2*τ + m*g*y + m*(u^2 + v^2))\n","end\n","M = Diagonal([1,1,1,1,0])\n","f = ODEFunction(pendulum, mass_matrix=M)\n","problem = ODEProblem(f,u₀,tspan,(ℓ,1,1))\n","solution = solve(problem, Rodas5(), reltol=1e-8, abstol=1e-8)\n","plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using ModelingToolkit, DifferentialEquations\n","θ₀ = π/3; ℓ = 1; tspan = (0,30.0)\n","u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, 0]\n","M = Diagonal([1,1,1,1,0])\n","function pendulum(dU,U,p,t)\n"," x,y,u,v,τ = U; ℓ,g,m = p\n"," dU .= (u, v, -τ*x/m, -τ*y/m + g, x^2 + y^2 - ℓ^2)\n","end\n","f = ODEFunction(pendulum, mass_matrix=M)\n","problem = ODEProblem(f, u₀, tspan, [ℓ,1,1])\n","sys = modelingtoolkitize(problem)\n","pendulum_sys = structural_simplify(dae_index_lowering(sys))\n","problem = ODAEProblem(pendulum_sys, [], tspan)\n","solution = solve(problem, Tsit5(), abstol=1e-8, reltol=1e-8)\n","plot(solution, idxs=(1,3), yflip=true, aspect_ratio=:equal)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Gauss–Legendre–Lobatto orthogonal collocation** The following code solves the pendulum problem using orthogonal collocation."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FastGaussQuadrature\n","function differentiation_matrix(n,Δt=1) \n"," nodes, _ = gausslobatto(n+1)\n"," t = (nodes[2:end].+1)/2\n"," A = t.^(0:n-1)'.*(1:n)'\n"," B = t.^(1:n)'\n"," (A/B)/Δt, t\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function pendulumGLL(U,(U₀,D,n))\n"," x,y,u,v,τ = U[1:n],U[n+1:2n],U[2n+1:3n],U[3n+1:4n],U[4n+1:5n]\n"," x₀,y₀,u₀,v₀,_ = U₀\n"," ℓ,g,m = (1,1,1)\n"," [D(x,x₀) .- u;\n"," D(y,y₀) .- v;\n"," D(u,u₀) .+ τ.*x/m;\n"," D(v,v₀) .+ τ.*y/m .- g;\n"," x.^2 + y.^2 .- ℓ^2]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using NonlinearSolve\n","θ₀ = π/3; ℓ = 1; tspan = (0,30.0)\n","u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, 0]\n","n = 5; N = 100; Δt = 30/N\n","M,_ = differentiation_matrix(n,Δt) \n","D = (u,u₀) -> M*(u .- u₀)\n","u = u₀\n","for i = 2:N\n"," problem = NonlinearProblem(pendulumGLL,[ones(n)*u₀'...],(u₀,D,n))\n"," solution = solve(problem,NewtonRaphson(),abstol=1e-12)\n"," u = [u reshape(solution.u,5,n)'] \n"," u₀ = u[:,end]\n","end\n","plot(u[1,:], u[2,:], yflip=true, aspect_ratio=:equal)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 13: Parabolic Equations\n**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Δx = .01; Δt = .01; L = 2; ν = Δt/Δx^2; uₗ = 0; uᵣ = 0;\n","x = -L:Δx:L; m = length(x) \n","u = (abs.(x).<1)\n","u[1] += 2ν*uₗ; u[m] += 2ν*uᵣ \n","D = Tridiagonal(ones(m-1), -2ones(m), ones(m-1))\n","D[1,2] = 0; D[m,m-1] = 0\n","A = I - ν*D \n","for i in 1:20\n"," u = A\\u\n","end\n","plot(x,u)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Δx = .01; Δt = .01; L = 2; ν = Δt/Δx^2\n","x = -L:Δx:L; m = length(x) \n","u = float.(abs.(x).<1)\n","D = Tridiagonal(ones(m-1), -2ones(m), ones(m-1))\n","D[1,2] = 2; D[m,m-1] = 2\n","A = 2I + ν*D \n","B = 2I - ν*D \n","anim = @animate for i in 1:40\n"," plot(x,u, legend=:none, ylims=(0,1))\n"," u .= B\\(A*u)\n","end\n","gif(anim, \"heat_cn.gif\", fps = 5)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Porous medium equation.** We'll now solve the porous medium equation $u_t = (u^2u_x)_x$ using an adaptive-step BDF routine."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Sundials\n","m = 400; L = 2; x = LinRange(-L,L,m); Δx = x[2]-x[1]\n","α = (u -> u.^2)\n","Du(u,Δx,t) = [0;diff(α((u[1:end-1]+u[2:end])/2).*diff(u))/Δx^2;0]\n","u₀ = (abs.(x).<1)\n","problem = ODEProblem(Du,u₀,(0,2),Δx)\n","method = CVODE_BDF(linear_solver=:Band, jac_lower=1, jac_upper=1)\n","solution = solve(problem, method);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["anim = @animate for t in LinRange(0,2,200)\n"," plot(x,solution(t),legend=:none,fill=(0,0.4,:red),ylims=(0,1))\n","end\n","gif(anim, \"porous.gif\", fps = 15)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Interact\n","@manipulate for t∈slider(0:0.01:2; value=0, label=\"time\")\n"," plot(x,solution(t), fill = (0, 0.4, :red))\n"," plot!(ylims=(0,1),legend=:none)\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Chapter 16: Fourier Spectral Methods\n**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].$$"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW\n","m = 256; ℓ = 4\n","ξ² = (fftshift(-(m-1)÷2:m÷2)*(2π/ℓ)).^2\n","u = (t,u₀) -> real(ifft(exp.(-ξ²*t).*fft(u₀)))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["x = LinRange(-ℓ/2,ℓ/2,m+1)[2:end]\n","u₀ = (abs.(x).<1)\n","plot(x,u(0.1,u₀))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["F = plan_fft(u₀)\n","u = (t,u₀) -> real(F\\(exp.(-ξ²*t).*(F*u₀)))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Kuramoto–Sivashinsky equation.** Let's solve the KSE using an IMEX method."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["#= In moving from Julia 1.7 to 1.9 this code no longer works.\n","using FFTW, DifferentialEquations, SciMLOperators\n","ℓ = 100; T = 200.0; n = 512\n","x = LinRange(-ℓ/2,ℓ/2,n+1)[2:end]\n","u₀ = exp.(-x.^2)\n","iξ = im*[0:(n÷2);(-n÷2+1):-1]*2π/ℓ\n","F = plan_fft(u₀)\n","L = DiagonalOperator(-iξ.^2-iξ.^4)\n","N = (u,_,_) -> -0.5iξ.*(F*((F\\u).^2))\n","problem = SplitODEProblem(L,N,F*u₀,(0,T))\n","method = KenCarp3(linsolve=KrylovJL(),autodiff=false)\n","solution = solve(problem,method,reltol=1e-12)\n","u = t -> real(F\\solution(t))\n","=#"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["#=\n","using ImageShow,ColorSchemes\n","s = hcat(u.(LinRange(0,T,500))...)\n","get(colorschemes[:binary],abs.(s),:extrema)\n","=#"]},{"cell_type": "markdown", "metadata": {}, "source": ["**Navier–Stokes equation.** We solve the Navier–Stokes equation in three parts: define the functions, initialize the variables, and iterate over time."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Δᵒ(Q,step=1) = Q - circshift(Q,(step,0))\n","flux(Q,c) = c.*Δᵒ(Q) - 0.5c.*(1 .- c).*(Δᵒ(Q,1)+Δᵒ(Q,-1))\n","H = (u,v,iξ₁,iξ₂) -> F*((F\\u).*(F\\(iξ₁.*u)) + (F\\v).*(F\\(iξ₂.*u)))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW\n","ℓ = 2; n = 128; ϵ = 0.001; Δt = 0.001; Δx = ℓ/n\n","x = LinRange(Δx,ℓ,n)'; y = x'\n","Q = (@. 0.5(1 + tanh(10 - 20abs(1 - 2y/ℓ)))).*ones(1,n)\n","uⁿ = Q .* (1 .+ 0.5sin.(ℓ*π*x))\n","vⁿ = zeros(n,n) \n","F = plan_fft(uⁿ)\n","uⁿ, vⁿ = F*uⁿ, F*vⁿ \n","uᵒ, vᵒ = uⁿ, vⁿ\n","iξ₁ = im*fftfreq(n,n)'*(2π/ℓ)\n","iξ₂ = transpose(iξ₁)\n","ξ² = iξ₁.^2 .+ iξ₂.^2\n","H₁ⁿ, H₂ⁿ = H(uⁿ,vⁿ,iξ₁,iξ₂), H(vⁿ,uⁿ,iξ₂,iξ₁)\n","M₊ = (1/Δt .+ (ϵ/2)*ξ²)\n","M₋ = (1/Δt .- (ϵ/2)*ξ²);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["for i = 1:1200\n"," Q -= flux(Q,(Δt/Δx).*real(F\\vⁿ)) + flux(Q',(Δt/Δx).*real(F\\uⁿ)')'\n"," H₁ⁿ⁻¹, H₂ⁿ⁻¹ = H₁ⁿ, H₂ⁿ \n"," H₁ⁿ, H₂ⁿ = H(uⁿ,vⁿ,iξ₁,iξ₂), H(vⁿ,uⁿ,iξ₂,iξ₁)\n"," uᵒ = uⁿ - uᵒ + (-1.5H₁ⁿ + 0.5H₁ⁿ⁻¹ + M₊.*uⁿ)./M₋\n"," vᵒ = vⁿ - vᵒ + (-1.5H₂ⁿ + 0.5H₂ⁿ⁻¹ + M₊.*vⁿ)./M₋\n"," ϕ = (iξ₁.*uᵒ + iξ₂.*vᵒ)./(ξ² .+ (ξ².≈ 0)) \n"," uⁿ, vⁿ = uᵒ - iξ₁.*ϕ, vᵒ - iξ₂.*ϕ\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using ColorSchemes\n","imresize(get(ColorSchemes.acton, clamp01.(Q)),ratio=3)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n# Part 4: Solutions\n\n## Numerical Linear Algebra\n**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](http://oeis.org/A055165).) 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$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using LinearAlgebra, Plots\n","N = 10000\n","n = [sum([!(det(rand(Bool,d,d))≈0) \n"," for i in 1:N]) for d in 1:20]\n","plot(n/N,marker=:o)\n","p = [1,6,174,22560,12514320,28836612000,270345669985440,10160459763342013440]\n","r = @. 2^(log2(p)-(1:8)^2)\n","plot!(r,marker=:square,opacity=0.5)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = {{0,1,0,1},{1,0,2,1},{0,2,0,0},{1,1,0,0}};\n","W = Inverse[IdentityMatrix[4] - z*A];\n","Simplify[W[[1,3]]]\n","Series[%,{z,0,10}]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**2.4. Good Will Hunting problem.** The following function computes the number of walks for $n$-step walks for the graph
\"good"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = [0 1 0 1;1 0 2 1;0 2 0 0;1 1 0 0]\n","walks(A,i,j,N) = [(A^n)[i,j] for n in 1:N]\n","walks(A,1,3,9)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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}$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function detx(A)\n"," L,U,P = lu(A)\n"," s = 1\n"," for i in 1:length(P)\n"," m = findfirst(P.==i)\n"," i!=m && ( s *= -1; P[[m,i]]=P[[i,m]] )\n"," end\n"," s * prod(diag(U))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = rand(20,20)\n","detx(A) ≈ det(A)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function rcuthillmckee(A)\n"," r = sortperm(vec(sum(A.!=0,dims=2)))\n"," p = Int64[]\n"," while ~isempty(r)\n"," q = [popfirst!(r)]\n"," while ~isempty(q)\n"," q₁ = popfirst!(q)\n"," append!(p,q₁) \n"," k = findall(!iszero, A[q₁,r])\n"," append!(q,splice!(r,k))\n"," end\n"," end\n"," reverse(p)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using SparseArrays\n","A = sparse(Symmetric(sprand(1000,1000,0.003)))\n","p = rcuthillmckee(A)\n","plot(spy(A), spy(A[p,p]), colorbar = false)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**2.6. Dolphins of Doubtful Sound.** We'll reuse the code [above](#dolphins_graph) used to draw the original graph of the dolphins."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DelimitedFiles, SparseArrays, Graphs, GraphPlot\n","ij = readdlm(download(bucket*\"dolphins.csv\"),',',Int)\n","A = sparse(ij[:,1],ij[:,2],one.(ij[:,1]))\n","gplot(Graph(A),layout=circular_layout)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["p = rcuthillmckee(A)\n","gplot(Graph(A[p,p]),layout=circular_layout)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**2.8. Stigler diet problem.** Let's solve the Stigler diet problem. We'll use the naïve [`simplex`](#simplex) function defined above."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using CSV, DataFrames\n","diet = DataFrame(CSV.File(download(bucket*\"diet.csv\")))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = Array(diet[2:end,4:end])'\n","c = ones.(size(A,2))\n","b = Array(diet[1,4:end])\n","food = diet[2:end,1]\n","solution = simplex(b,A',c)\n","print(\"foods: \", food[solution.y .!= 0], \"\\n\")\n","print(\"daily cost: \", solution.z)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using JuMP, GLPK\n","model = Model(GLPK.Optimizer)\n","@variable(model, x[1:size(A,2)] ≥ 0)\n","@objective(model, Min, c' * x)\n","@constraint(model, A * x .≥ b)\n","optimize!(model)\n","print(\"foods: \", food[JuMP.value.(x) .!= 0], \"\\n\")\n","print(\"daily cost: \", objective_value(model))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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}.$$"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function findpath(A,a,b)\n"," r = collect(1:size(A,2))\n"," q = [a]; p = [-1]; i = 0\n"," splice!(r,a)\n"," while i findfirst(actors.==x)[1]\n","actormovie[findpath(A, index(\"Emma Watson\"), index(\"Kevin Bacon\"))]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["actormovie[findpath(A, index(\"Bruce Lee\"), index(\"Tommy Wiseau\"))]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**2.10. Sparse matrix storage efficiency.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["d = 0.5; A = sprand(800,600,d)\n","Base.summarysize(A)/Base.summarysize(Matrix(A))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["X = load(download(bucket*\"laura.png\")) .|> Gray\n","m,n = size(X)\n","blur = x -> exp(-(x/20).^2/2)\n","A = [blur(i-j) for i=1:m, j=1:m]; A ./= sum(A,dims=2)\n","B = [blur(i-j) for i=1:n, j=1:n]; B ./= sum(B,dims=1)\n","N = 0.01*rand(m,n)\n","Y = A*X*B + N;"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["α = 0.05\n","X₁ = A\\Y/B\n","X₂ = (A'*A+α^2*I)\\A'*Y*B'/(B*B'+α^2*I)\n","X₃ = pinv(A,α)*Y*pinv(B,α)\n","Gray.([X Y X₁ X₂ X₃])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["vandermonde(t,n) = vec(t).^(0:n-1)'\n","build_poly(c,X) = vandermonde(X,length(c))*c"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function solve_filip(x,y,n)\n"," V = vandermonde(x, n)\n"," c = Array{Float64}(undef, 3, n)\n"," c[1,:] = (V'*V)\\(V'*y)\n"," Q,R = qr(V)\n"," c[2,:] = R\\(Matrix(Q)'*y)\n"," c[3,:] = pinv(V,1e-9)*y\n"," r = [norm(V*c[i,:]-y) for i in 1:3]\n"," return(c,r)\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["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?"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DelimitedFiles\n","data = readdlm(download(bucket*\"filip.csv\"),',',Float64)\n","coef = readdlm(download(bucket*\"filip-coeffs.csv\"),',')\n","(x,y) = (data[:,2],data[:,1])\n","β,r = solve_filip(x, y, 11)\n","X = LinRange(minimum(x),maximum(x),200)\n","Y = [build_poly(β[i,:],X) for i in 1:3]\n","plot(X,Y); scatter!(x,y,opacity=0.5)\n","[coef β']"]},{"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Statistics\n","zscore(X,x=X) = (X .- mean(x))/std(x)\n","c,r = solve_filip(zscore(x), zscore(y), 11)\n","Y = [build_poly(c[i,:],zscore(X,x))*std(y).+mean(y) for i in 1:3]\n","plot(X,Y); scatter!(x,y,opacity=0.5)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using CSV, DataFrames, Dates\n","data = download(bucket*\"dailytemps.csv\")\n","df = DataFrame(CSV.File(data))\n","day = Dates.value.(df[:,:date] .- df[1,:date])/365\n","u = df[:,:temperature]\n","tempsmodel(t) = [sin.(2π*t) cos.(2π*t) one.(t)]\n","c = tempsmodel(day)\\u\n","scatter(df[:,:date],u,alpha=0.3)\n","plot!(df[:,:date],tempsmodel(day)*c,width=3)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**3.8. Image recognition.** We'll practice identifying handwritten digits using the MNIST database. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using MLDatasets, Arpack, Images\n","image_train, label_train = MNIST(split=:train)[:]\n","image_train = reshape(permutedims(image_train,[3,2,1]),60000,:)\n","V = [( D = image_train[label_train.==i,:];\n"," svds(D,nsv=12)[1].V ) for i ∈ 0:9];"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["pix = -abs.(reshape(V[3+1],28,:))\n","rescale = scaleminmax(minimum(pix), maximum(pix))\n","pix .|> rescale .|> Gray"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["image_test, label_test = MNIST(split=:test)[:]\n","image_test = reshape(permutedims(image_test,[3,2,1]),10000,:)\n","r = [( q = image_test*V[i]*V[i]' .- image_test;\n"," sum(q.^2,dims=2) ) for i∈1:10]\n","r = hcat(r...) \n","prediction = [argmin(r[i,:]) for i∈1:10000] .- 1\n","[sum(prediction[label_test.== i].== j) for j∈0:9, i∈0:9]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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`](#get_names) and [`get_adjacency_matrix`](#get_adjacency_matrix) developed earlier."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["actors = get_names(\"actors\")\n","genres = get_names(\"genres\")\n","A = get_adjacency_matrix(\"movie-genre\")\n","A = (A*Diagonal(1 ./ sum(A,dims=1)[:]))\n","B = get_adjacency_matrix(\"actor-movie\");"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Arpack\n","X,_ = svds(A*B,nsv=10)\n","U = X.U; Σ = Diagonal(X.S); Vᵀ = X.Vt\n","Q = Vᵀ./sqrt.(sum(Vᵀ.^2,dims=1));"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["index = x -> findfirst(actors.==x)[1]\n","q = Q[:,index(\"Steve Martin\")][:]\n","z = Q'*q\n","r = sortperm(-z)\n","actors[r[1:10]]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["p = U*Σ*q\n","r = sortperm(-p)\n","[genres[r] p[r]/sum(p)]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**3.10. Multilateration.** We use ordinary least squares and total least squares to solve a multilateration problem. The function [`tls`](#tls) is defined above."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["X = [3 3 12; 1 15 14; 10 2 13; 12 15 14; 0 11 12]\n","reference = X[1,:]; X .-= reference'\n","A = [2 2 -2].*X; b = (X.^2)*[1; 1; -1]\n","xₒₗₛ = A\\b + reference\n","xₜₗₛ = tls(A,b) + reference\n","xₒₗₛ, xₜₗₛ"]},{"cell_type": "markdown", "metadata": {}, "source": ["**3.11. ShotSpotter.** Let's modify the previous solution for this multilateration problem. We'll first download the data set."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using CSV, DataFrames\n","df = DataFrame(CSV.File(download(bucket*\"shotspotter.csv\")))\n","X = Array(df[1:end-1,:]); xₒ = Array(df[end,:]);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["c² = 328.9^2\n","reference = X[1,:]; X .-= reference'\n","A = [2 2 2 -2c²].*X; b = (X.^2)*[1; 1; 1; -c²]\n","xₒₗₛ = A\\b + reference\n","error = xₒₗₛ - xₒ"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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?"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 8\n","E = collect(Iterators.flatten([eigvals(randn(n,n)) for i∈1:2000]))\n","scatter(E,mc=:black,ma=0.05,legend=nothing,aspect_ratio=:equal)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function rayleigh(A)\n"," x = randn(size(A,1),1)\n"," while true\n"," x = x/norm(x) \n"," ρ = (x'*A*x)[1]\n"," M = A - ρ*I\n"," abs(cond(M, 1)) < 1e12 ? x = M\\x : return(ρ,x)\n"," end\n","end\n","A = [2 3 6 4; 3 0 3 1; 6 3 8 8; 4 1 8 2]\n","rayleigh(A)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function implicitqr(A)\n"," n = size(A,1)\n"," tolerance = 1E-12\n"," H = Matrix(hessenberg(A).H)\n"," while true\n"," if abs(H[n,n-1]) sortperm ][1:10]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function randomizedsvd(A,k)\n"," Z = rand(size(A,2),k);\n"," Q = Matrix(qr(A*Z).Q)\n"," for i in 1:3\n"," Q = Matrix(qr(A'*Q).Q)\n"," Q = Matrix(qr(A*Q).Q)\n"," end\n"," W = svd(Q'*A)\n"," return((Q*W.U,W.S,W.V)) \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Images\n","img = load(download(bucket*\"red-fox.jpg\"))\n","A = Float64.(Gray.(img))\n","U,S,V = randomizedsvd(A,10)\n","Gray.([A U*Diagonal(S)*V'])"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Arpack\n","@time randomizedsvd(A,10)\n","@time svds(A, nsv=10)\n","@time svd(A);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϵ = 1 .- sqrt.(cumsum(S.^2))/norm(A)\n","plot(ϵ,marker=:circle)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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$?\""]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Arpack\n","m = 5000\n","A = [1/(1 + (i+j-1)*(i+j)/2 - j) for i∈1:m, j∈1:m]\n","svds(A,nsv=1)[1].S[1]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using SparseArrays, LinearAlgebra\n","⊗(x,y) = kron(x,y); ϕ = x -> x-x^2\n","n = 50 ; x = (1:n)/(n+1); Δx = 1/(n+1)\n","J = sparse(I, n, n)\n","D = spdiagm(-1 => ones(n-1), 0 => -2ones(n), 1 => ones(n-1) )\n","A = ( D⊗J⊗J + J⊗D⊗J + J⊗J⊗D ) / Δx^2\n","f = [ϕ(x)*ϕ(y) + ϕ(x)*ϕ(z) + ϕ(y)*ϕ(z) for x∈x, y∈x, z∈x][:]\n","uₑ = [ϕ(x)*ϕ(y)*ϕ(z)/2 for x∈x, y∈x, z∈x][:];"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function stationary(A,b,ω=0,n=400)\n"," ϵ = []; u = zero.(b)\n"," P = (1-ω)*sparse(Diagonal(A)) + ω*sparse(LowerTriangular(A))\n"," for i=1:n\n"," u += P\\(b-A*u)\n"," append!(ϵ,norm(u - uₑ,1))\n"," end\n"," return(ϵ)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function conjugategradient(A,b,n=400)\n"," ϵ = []; u = zero.(b) \n"," p = r = b-A*u\n"," for i=1:n\n"," Ap = A*p\n"," α = (r'*p)/(Ap'*p)\n"," u += α.*p; r -= α.*Ap\n"," β = (r'*Ap)/(Ap'*p)\n"," p = r - β.*p\n"," append!(ϵ,norm(u - uₑ,1))\n"," end\n"," return(ϵ)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϵ = zeros(400,4)\n","@time ϵ[:,1] = stationary(A,-f,0)\n","@time ϵ[:,2] = stationary(A,-f,1)\n","@time ϵ[:,3] = stationary(A,-f,1.9)\n","@time ϵ[:,4] = conjugategradient(A,-f)\n","method = [\"Jacobi\" \"Gauss-Seidel\" \"SOR\" \"Conjugate Gradient\"]\n","plot(ϵ,yaxis=:log,labels=method,bglegend=RGBA(1,1,1,0.7))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["k = 1:120; ([one.(k) k]\\log10.(ϵ[k,:]))[2,:]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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}$?\""]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Primes, LinearAlgebra, SparseArrays, IterativeSolvers\n","n = 20000\n","d = 2 .^ (0:14); d = [-d;d]\n","P = Diagonal(primes(224737))\n","B = [ones(n - abs(d)) for d∈d]\n","A = sparse(P) + spdiagm(n ,n, (d .=> B)...)\n","b = zeros(n); b[1] = 1\n","cg(A, b; Pl=P)[1]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**6.1. Radix-3 FFT.** The radix-3 FFT is similar to the [radix-2 FFT](#radix2fft). We'll verify that the code is correct by comparing it with a built-in FFT"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function fftx3(c)\n"," n = length(c)\n"," ω = exp(-2im*π/n)\n"," if mod(n,3) == 0 \n"," k = collect(0:n/3-1)\n"," u = [transpose(fftx3(c[1:3:n-2]));\n"," transpose((ω.^k).*fftx3(c[2:3:n-1]));\n"," transpose((ω.^2k).*fftx3(c[3:3:n]))] \n"," F = exp(-2im*π/3).^([0;1;2]*[0;1;2]')\n"," return(reshape(transpose(F*u),:,1))\n"," else \n"," F = ω.^(collect(0:n-1)*collect(0:n-1)')\n"," return(F*c)\n"," end \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW\n","v = rand(24,1)\n","[fft(v) fftx3(v)]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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](https://en.wikipedia.org/wiki/RSA_numbers#RSA-129). "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW\n","function multiply(p_,q_)\n"," p = [parse.(Int,split(reverse(p_),\"\"));zeros(length(q_),1)]\n"," q = [parse.(Int,split(reverse(q_),\"\"));zeros(length(p_),1)]\n"," pq = Int.(round.(real.(ifft(fft(p).*fft(q)))))\n"," carry = pq .÷ 10\n"," while any(carry.>0)\n"," pq -= carry*10 - [0;carry[1:end-1]]\n"," carry = pq .÷ 10\n"," end\n"," n = findlast(x->x>0, pq)\n"," return(reverse(join(pq[1:n[1]])))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["p = \"32769132993266709549961988190834461413177642967992942539798288533\"\n","q = \"3490529510847650949147849619903898133417764638493387843990820577\"\n","multiply(p,q), parse(BigInt, p)*parse(BigInt, q)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Random\n","p = randstring('0':'9', 100000)\n","q = randstring('0':'9', 100000)\n","@time multiply(p,q)\n","@time parse(BigInt, p)*parse(BigInt, q);"]},{"cell_type": "markdown", "metadata": {}, "source": ["**6.3. Fast discrete cosine transform.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function dctII(f)\n"," n = size(f,1)\n"," ω = exp.(-0.5im*π*(0:n-1)/n)\n"," return real(ω.*fft(f[[1:2:n; n-mod(n,2):-2:2],:],1))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function idctII(f)\n"," n = size(f,1)\n"," ω = exp.(-0.5im*π*(0:n-1)/n)\n"," f[1,:] = f[1,:]/2\n"," f[[1:2:n; n-mod(n,2):-2:2],:] = 2*real(ifft(f./ω,1))\n"," return f\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["dct2(f) = dctII(dctII(f')')\n","idct2(f) = idctII(idctII(f')')"]},{"cell_type": "markdown", "metadata": {}, "source": ["**6.4. DCT image compression.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function dctcompress2(A,d)\n"," B = dct(Float64.(A))\n"," m,n = size(A) \n"," m₀,n₀ = sqrt(d).*(m,n) .|> floor .|> Int\n"," B₀ = B[1:m₀,1:n₀]\n"," A₀ = idct([B₀ zeros(m₀,n-n₀); zeros(m-m₀,n)])\n"," A₀ = A₀ |> clamp01! .|> Gray\n"," sizeof(B₀)/sizeof(B), sqrt(1-(norm(B₀)/norm(B))^2), A₀\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = Gray.(load(download(bucket*\"laura_square.png\")))\n","[A dctcompress2(A,0.05)[3]]"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Numerical Analysis\n**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["aitken₁(x₁,x₂,x₃) = x₃ - (x₃-x₂)^2/(x₃ - 2x₂ + x₁)\n","aitken₂(x₁,x₂,x₃) = (x₁*x₃ - x₂^2)/(x₃ - 2x₂ + x₁)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 60000\n","p = cumsum([(-1)^i*4/(2i+1) for i=0:n])\n","p₁ = aitken₁.(p[1:n-2],p[2:n-1],p[3:n])\n","p₂ = aitken₂.(p[1:n-2],p[2:n-1],p[3:n])\n","plot(abs.(π.-p)); plot!(abs.(π.-p₂)); plot!(abs.(π.-p₁))\n","plot!(xaxis=:log,yaxis=:log)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations\n","function homotopy(f,df,x)\n"," dxdt(x,p,t) = -df(x)\\p\n"," sol = solve(ODEProblem(dxdt,x,(0.0,1.0),f(x)))\n"," sol.u[end]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function newton(f,df,x)\n"," for i in 1:100\n"," Δx = -df(x)\\f(x)\n"," norm(Δx) > 1e-8 ? x += Δx : return(x)\n"," end\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["f = z -> ((x,y)=tuple(z...);\n"," [(x^2+y^2)^2-2(x^2-y^2); (x^2+y^2-1)^3-x^2*y^3])\n","df = z -> ((x,y)=tuple(z...);\n"," [4x*(x^2+y^2-1) 4y*(x^2+y^2+1);\n"," 6x*(x^2+y^2-1)^2-2x*y^3 6y*(x^2+y^2-1)^2-3x^2*y^2])"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["display(homotopy(f,df,[1,1]))\n","display(newton(f,df,[1,1]))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function ⊕(P,Q)\n"," a = 0\n"," r = BigInt(2)^256 - 2^32 - 977\n"," if P[1]==Q[1]\n"," d = invmod(2*P[2], r)\n"," λ = mod((3*powermod(P[1],2,r)+a)*d, r)\n"," else\n"," d = invmod(Q[1]-P[1], r)\n"," λ = mod((Q[2]-P[2])*d, r)\n"," end\n"," x = mod(powermod(λ, 2, r) - P[1] - Q[1], r)\n"," y = mod(-λ*(x-P[1])-P[2], r)\n"," [x;y]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Base.isodd(m::BigInt) = ((m&1)==1)\n","function dbl_add(m,P)\n"," if m > 1\n"," Q = dbl_add(m>>1,P)\n"," return isodd(m) ? (Q⊕Q)⊕P : Q⊕Q\n"," else\n"," return P\n"," end\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["P₁ = big\"0x79BE667EF9DCBBAC55A06295CE87\n"," 0B07029BFCDB2DCE28D959F2815B16F81798\"\n","P₂ = big\"0x483ADA7726A3C4655DA4FBFC0E11\n"," 08A8FD17B448A68554199C47D08FFB10D4B8\"\n","P = [P₁; P₂]\n","m, n = 1829442, 3727472\n","mP = dbl_add(m,P)\n","nmP = dbl_add(n,mP)\n","nP = dbl_add(n,P)\n","print(\"Alice's private key:\\n$m\\n\\n\")\n","print(\"Bob's private key:\\n$n\\n\\n\")\n","print(\"Alice's public key:\\n$mP\\n\\n\")\n","print(\"Bob's public key:\\n$nP\\n\\n\")\n","print(\"Alice and Bob's shared secret key:\\n$nmP\")"]},{"cell_type": "markdown", "metadata": {}, "source": ["**9.2. Periodic parametric splines.** We modify the code [`spline_natural`](#spline_natural) (above) to make a generate a spine with periodic boundary conditions. The function `evaluate_spline` is duplicated from the code above."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function spline_periodic(x,y)\n"," h = diff(x) \n"," d = 6*diff(diff([y[end-1];y])./[h[end];h])\n"," α = h[1:end-1]\n"," β = 2*(h+circshift(h,1))\n"," C = Matrix(SymTridiagonal(β,α))\n"," C[1,end]=h[end]; C[end,1]=h[end] \n"," m = C\\d\n"," return([m;m[1]])\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function evaluate_spline(xᵢ,yᵢ,m,n)\n"," h = diff(xᵢ)\n"," B = yᵢ[1:end-1] .- m[1:end-1].*h.^2/6\n"," A = diff(yᵢ)./h - h./6 .*diff(m)\n"," x = range(minimum(xᵢ),maximum(xᵢ),length=n+1)\n"," i = sum(xᵢ' .≤ x,dims=2)\n"," i[end] = length(xᵢ)-1 \n"," y = @. (m[i]*(xᵢ[i+1]-x)^3 + m[i+1]*(x-xᵢ[i])^3)/(6h[i]) +\n"," A[i]*(x-xᵢ[i]) + B[i]\n"," return(x,y)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 5; nx = 30\n","x = rand(n); y = rand(n)\n","x = [x;x[1]]; y = [y;y[1]]\n","t = [0;cumsum(sqrt.(diff(x).^2 + diff(y).^2))]\n","X = evaluate_spline(t,x,spline_periodic(t,x),nx*n)\n","Y = evaluate_spline(t,y,spline_periodic(t,y),nx*n)\n","scatter(x,y); plot!(X[2],Y[2],legend=false)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 20; N = 200\n","x = LinRange(-1,1,n)\n","y = float(x .> 0);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϕ₁(x,a) = abs.(x.-a).^3\n","ϕ₂(x,a) = exp.(-20(x.-a).^2)\n","ϕ₃(x,a) = x.^a"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["X = LinRange(-1,1,N)\n","interp(ϕ,a) = ϕ(X,a')*(ϕ(x,a')\\y)\n","Y₁ = interp(ϕ₁,x)\n","Y₂ = interp(ϕ₂,x)\n","Y₃ = interp(ϕ₃,(0:n-1))\n","scatter(x,y,seriestype = :scatter, marker = :o, legend = :none)\n","plot!(X,[Y₁,Y₂,Y₃],ylim=[-0.5,1.5])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function collocation_solve(L,f,bc,x)\n"," h = x[2]-x[1]\n"," S = L(x)*([1 -1/2 1/6; -2 0 2/3; 1 1/2 1/6]./[h^2 h 1])'\n"," d = [bc[1]; f(x); bc[2]]\n"," A = Matrix(Tridiagonal([S[:,1];0], [0;S[:,2];0], [0;S[:,3]]))\n"," A[1,1:3], A[end,end-2:end] = [1 4 1]/6, [1 4 1]/6 \n"," return(A\\d)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function collocation_build(c,x,N)\n"," X = LinRange(x[1],x[end],N)\n"," h = x[2] - x[1]\n"," i = Int32.(X .÷ h .+ 1); i[N] = i[N-1]\n"," C = [c[i] c[i.+1] c[i.+2] c[i.+3]]'\n"," B = (x->[(1-x).^3;4-3*(2-x).*x.^2;4-3*(1+x).*(1-x).^2;x.^3]/6)\n"," Y = sum(C.*hcat(B.((X.-x[i])/h)...),dims=1)\n"," return(Array(X),reshape(Y,(:,1)))\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["Now, we can solve the Bessel equation $xu''+u'+xu =0$ with boundary conditions $u(0)=1$ and $u(b)=0$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Roots, SpecialFunctions\n","n = 20; N = 141\n","L = (x -> [x one.(x) x])\n","f = (x -> zero.(x) )\n","b = fzero(besselj0, 11)\n","x = range(0,b,length=n)\n","c = collocation_solve(L,f,[1,0],x)\n","X,Y = collocation_build(c,x,70)\n","plot(X,[Y besselj0.(X)])"]},{"cell_type": "markdown", "metadata": {}, "source": ["Finally, we'll explore the error and convergence rate."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["N = 10*2 .^(1:7); ϵ = []\n","for n in N\n"," x = LinRange(0,b,n) \n"," c = collocation_solve(L,f,[1,0],x)\n"," X,Y = collocation_build(c,x,n)\n"," append!(ϵ, norm(Y-besselj0.(X)) / n)\n","end\n","slope = ([log.(N) one.(N)]\\log.(ϵ))[1]\n","s = \"slope = \"*string.(round.(slope,digits=2))\n","plot(N, ϵ, xaxis=:log, yaxis=:log, marker=:o, label=s)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**10.3. Fractional derivatives.** We'll plot the fractional derivatives for a function."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW\n","n = 256; ℓ = 2\n","x = (0:n-1)/n*ℓ .- ℓ/2\n","ξ = [0:(n/2-1); (-n/2):-1]*(2π/ℓ)\n","f₁ = exp.(-16*x.^2)\n","f₂ = sin.(π*x)\n","f₃ = x.*(1 .- abs.(x))\n","deriv(f,p) = real(ifft((im*ξ).^p.*fft(f)))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Interact\n","func = Dict(\"Gaussian\"=>f₁,\"sine\"=>f₂,\"spline\"=>f₃)\n","@manipulate for f in togglebuttons(func; label=\"function\"), \n"," p in slider(0:0.01:2; value=0, label=\"derivative\")\n"," plot(x,deriv(f,p),legend=:none)\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**10.4. Handwriting classification.** We'll use Flux to train a convolutional neural net using MNIST data and then test the model."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using MLDatasets, Flux\n","model = Chain(\n"," Conv((5,5), 1=>6, tanh, pad=SamePad()),\n"," MeanPool((2,2)),\n"," Conv((5,5), 6=>16, tanh),\n"," MeanPool((2,2)),\n"," Conv((5,5), 16=>120, tanh),\n"," Flux.flatten,\n"," Dense(120, 84, tanh),\n"," Dense(84, 10))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["image_train, label_train = MLDatasets.MNIST(split=:train)[:]\n","image_test, label_test = MLDatasets.MNIST(split=:test)[:]\n","image_train = Flux.unsqueeze(image_train, 3) \n","image_test = Flux.unsqueeze(image_test, 3) \n","label_train = Flux.onehotbatch(label_train, 0:9)\n","label_test = Flux.onehotbatch(label_test, 0:9);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["loss(x,y) = Flux.Losses.logitcrossentropy(model(x),y)\n","parameters = Flux.params(model)\n","data = Flux.Data.DataLoader((image_train,label_train),batchsize=100)\n","optimizer = ADAM()"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using ProgressMeter\n","@showprogress for epochs = 1:5\n"," Flux.train!(loss, parameters, data, optimizer)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["accuracy(x,y) = sum(Flux.onecold(x) .== Flux.onecold(y))/size(y,2)\n","accuracy(model(image_test),label_test)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**10.5. Multilateration.** Let's modify the previous solution for this multilateration problem. We'll first download the data set."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using LsqFit\n","ϵ = (x,p) -> @. √((x[:,1]-p[1])^2 + (x[:,2]-p[2])^2) - (x[:,3]-p[3])\n","x = [3 3 12; 1 15 14; 10 2 13; 12 15 14; 0 11 12]\n","curve_fit(ϵ,x,zeros(5),zeros(3)).param"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using CSV, DataFrames\n","df = DataFrame(CSV.File(download(bucket*\"shotspotter.csv\")))\n","X = Array(df[1:end-1,:]); xₒ = Array(df[end,:]);"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using LsqFit\n","c = 328.6\n","ϵ = (x,p)->sqrt.(sum([(x[:,i].-p[i]).^2 for i∈1:3]))-c*(x[:,4].-p[4])\n","xₙₗₛ = curve_fit(ϵ,X,zero(X[:,1]),X[1,:]).param\n","error = xₙₗₛ - xₒ"]},{"cell_type": "markdown", "metadata": {}, "source": ["Let's also plot the solution. We'll first define a function to plot circles."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function plot_circle(x,y,r)\n"," t = LinRange(0,2π,100);\n"," plot!(x.+r*cos.(t), y.+r*sin.(t),color=\"black\",opacity=0.5);\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["r = 328.6*(X[:,end] .- xₙₗₛ[end])\n","scatter(X[:,1], X[:,2], color = \"black\")\n","[plot_circle(X[i,1],X[i,2],r[i]) for i in 1:length(r)]\n","scatter!([xₙₗₛ[1]], [xₙₗₛ[2]], color = \"red\")\n","plot!(legend=:none,aspect_ratio=:equal)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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)."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["d = [0,1,2,3]; n = length(d)\n","C = inv( d.^(0:n-1)' .// factorial.(0:n-1)' )\n","[C C*d.^n/factorial(n)]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**11.2. Richardson extrapolation.** The following code is an iterative version of the recursive [`richardson`](#richardson_extrapolation) function above:"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function richardson(f,x,m)\n"," D = []\n"," for i in 1:m\n"," append!(D, ϕ(f,x,2^i))\n"," for j in i-1:-1:1 \n"," D[j] = (4^(i-j)*D[j+1] - D[j])/(4^(i-j) - 1)\n"," end\n"," end\n"," D[1]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϕ = (f,x,n) -> (f(x+1/n) - f(x-1/n))/(2/n)\n","richardson(x->sin(x),0,4)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**11.3. Automatic differentiation.** Let's extend the [`Dual class`](#dualclass) 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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["struct Dual\n"," value\n"," deriv\n","end\n","\n","Dual(x) = Dual(x,1)\n","value(x) = x isa Dual ? x.value : x\n","deriv(x) = x isa Dual ? x.deriv : 0\n","\n","Base.:+(u, v) = Dual(value(u)+value(v), deriv(u)+deriv(v))\n","Base.:-(u, v) = Dual(value(u)-value(v), deriv(u)-deriv(v))\n","Base.:*(u, v) = Dual(value(u)*value(v), \n"," value(u)*deriv(v)+value(v)*deriv(u))\n","Base.:sin(u) = Dual(sin(value(u)), cos(value(u))*deriv(u))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Base.:sqrt(u) = Dual(sqrt(value(u)), deriv(u) / 2sqrt(value(u)))\n","Base.:/(u, v) = Dual(value(u)/value(v), \n"," (value(u)*deriv(v)-value(v)*deriv(u))/(value(v)*value(v)))\n","Base.:cos(u) = Dual(cos(value(u)), -sin(value(u))*deriv(u))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function get_zero(f,x)\n"," ϵ = 1e-12; δ = 1\n"," while abs(δ) > ϵ\n"," fx = f(Dual(x))\n"," δ = value(fx)/deriv(fx)\n"," x -= δ\n"," end\n"," return(x)\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["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}$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["get_zero(x->4sin(x)+sqrt(x),4)"]},{"cell_type": "markdown", "metadata": {}, "source": ["We can find a minimum or maximum of $4\\sin x + \\sqrt{x}$ by modifying Newton's method."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function get_extremum(f,x)\n"," ϵ = 1e-12; δ = 1\n"," while abs(δ)>ϵ\n"," fx = f(Dual(Dual(x)))\n"," δ = deriv(value(fx))/deriv(deriv(fx))\n"," x -= δ\n"," end\n"," return(x)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["get_extremum(x->4sin(x)+sqrt(x),4)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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}.$$"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function cauchyderivative(f, a, p, n = 20, ϵ = 0.1)\n"," ω = exp.(2π*im*(0:(n-1))/n)\n"," factorial(p)/(n*ϵ^p)*sum(@. f(a+ϵ*ω)/ω^p)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["f = x -> exp(x)/(cos(x)^3 + sin(x)^3)\n","cauchyderivative(f, 0, 6)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function gauss_legendre(n)\n"," x = -cos.((4*(1:n).-1)*π/(4n+2))\n"," Δ = one.(x)\n"," dPₙ = 0\n"," while(maximum(abs.(Δ))>1e-16)\n"," Pₙ, Pₙ₋₁ = x, one.(x)\n"," for k ∈ 2:n\n"," Pₙ, Pₙ₋₁ = ((2k - 1)*x.*Pₙ-(k-1)*Pₙ₋₁)/k, Pₙ \n"," end\n"," dPₙ = @. n*(x*Pₙ - Pₙ₋₁)/(x^2-1)\n"," Δ = Pₙ ./ dPₙ \n"," x -= Δ\n"," end\n"," return(x, @. 2/((1-x^2)*dPₙ^2))\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["f = (x-> 2sqrt(1-x^2))\n","x,w = gauss_legendre(10)\n","w ⋅ f.(x)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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.$$"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FastGaussQuadrature\n","ξ,w = gausshermite(40)\n","u₀ = x -> sin(x)\n","u = (t,x) -> w ⋅ u₀.(x.-2sqrt(t)*ξ)/sqrt(π)\n","x = LinRange(-12,12,100); plot(x,u.(1,x))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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})$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["mc_π(n) = sum(sum(rand(2,n).^2,dims=1).<1)/n*4"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m = 20; d = []; N = 2 .^ (1:20)\n","for n ∈ N\n"," append!(d,sum([abs(π - mc_π(n)) for i∈1:m])/m)\n","end\n","s = log.(N.^[0 1])\\log.(d)\n","scatter(N,d,xaxis=:log, yaxis=:log) \n","plot!(N,exp(s[1]).*N.^s[2])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**11.10 Orthogonal collocation** We'll solve $u'(t) = \\alpha u^2$ using the spectral method and pseudospectral method."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FastGaussQuadrature\n","function differentiation_matrix(n) \n"," nodes, _ = gausslobatto(n+1)\n"," t = (nodes[2:end].+1)/2\n"," A = t.^(0:n-1)'.*(1:n)'\n"," B = t.^(1:n)'\n"," A/B,t\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["n = 20\n","M,t = differentiation_matrix(n) \n","D = (u,u₀) -> M*(u .- u₀)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using NLsolve\n","u₀ = 1.0; α = 0.9\n","F = (u,u₀,α) -> D(u,u₀) .- α*u.^2\n","u = nlsolve(u->F(u,u₀,α),u₀*ones(n)).zero\n","plot([0;t], [u₀;u], marker=:circle, legend=false)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["N = 30; Δt = 1.0/N; n = 8\n","M,t = differentiation_matrix(n) \n","D = (u,u₀) -> M*(u .- u₀)/Δt\n","u₀ = 1.0; U = [u₀]; T = [0.0]; α = 0.9\n","for i = 0:N-1\n"," u = nlsolve(u->F(u,u₀,α),u₀*ones(n)).zero\n"," u₀ = u[end]\n"," append!(T,(i.+t)*Δt)\n"," append!(U,u)\n","end\n","plot(T, U, marker=:circle, legend=false)"]},{"cell_type": "markdown", "metadata": {}, "source": ["\n## Numerical Differential Equations\n**12.4. Runge–Kutta stability.** The following code plots the region of absolute stability for a Runge–Kutta method with tableau `A` and `b`."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["A = [0 0 0 0 0; \n"," 1/3 0 0 0 0; \n"," 1/6 1/6 0 0 0;\n"," 1/8 0 3/8 0 0;\n"," 1/2 0 -3/2 2 0];\n","b = [1/6 0 0 2/3 1/6];"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using LinearAlgebra, Plots\n","function rk_stability_plot(A,b)\n"," E = ones(length(b),1)\n"," r(λk) = abs.(1 .+ λk * b*((I - λk*A)\\E))\n"," x,y = LinRange(-4,4,100),LinRange(-4,4,100)\n"," s = reshape(vcat(r.(x'.+im*y)...),(100,100))\n"," contour(x,y,s,levels=[1],legend=false)\n","end\n","rk_stability_plot(A,b)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**12.7. Third-order IMEX coefficients.** We can determine the coefficients of a third-order IMEX method by inverting a system of linear equations."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["i = 0:3\n","a = ((-i)'.^i.//factorial.(i))\\[0;1;0;0]\n","b = ((-(i.+1))'.^i.//factorial.(i))\\[1;0;0;0]\n","[a b]"]},{"cell_type": "markdown", "metadata": {}, "source": ["**12.8. Predictor-corrector stability.** We'll use the [`multistepcoefficients`](#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$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function multistepcoefficients(m,n)\n"," s = length(m) + length(n) - 1\n"," A = (m.+1).^(0:s-1) .|> Rational\n"," B = (0:s-1).*((n.+1).^[0;0:s-2])\n"," c = -[A[:,2:end] B]\\ones(Int64,s)\n"," [1;c[1:length(m)-1]], c[length(m):end] \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Polynomials\n","function PECE(n,m)\n"," _,a = multistepcoefficients([0 1],hcat(1:n-1...))\n"," _,b = multistepcoefficients([0 1],hcat(0:n-1...))\n"," α(r) = a ⋅ r.^(1:n-1)/b[1]\n"," β(r) = b[2:end] ⋅ r.^(1:n-1)/b[1]\n"," z = [(c = [r-1; repeat([r + β(r)],m); α(r)];\n"," Polynomials.roots(Polynomial(c)))\n"," for r in exp.(im*LinRange(0,2π,200))]\n"," hcat(z/b[1]...)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["z = vcat([PECE(2,i)[:] for i in 0:4]...)\n","scatter(z,label=\"\",markersize=.5,aspect_ratio=:equal)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**12.9. Padé approximant.** We'll build a function to compute the coefficients of the Padé approximant to $log(r)$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function pade(a,m,n)\n"," A = Rational(1)*Matrix(I(m+n+1))\n"," A[:,m+2:end] = [i>j ? -a[i-j] : 0 for i∈1:m+n+1, j∈1:n]\n"," pq = A\\a[1:m+n+1]\n"," pq[1:m+1], [1; pq[m+2:end]]\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m = 3; n = 2\n","a = [0; ((-1).^(0:m+n).//(1:m+n+1))]\n","(p,q) = pade(a,m,n)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["S = n -> [(-1).^(i+j)*binomial(j,i) for i∈0:n, j∈0:n]\n","S(m)*p, S(n)*q"]},{"cell_type": "markdown", "metadata": {}, "source": ["**12.13. SIR solution.** Let's solve the susceptible-infected-recovered (SIR) model for infectious diseases using a general ODE solver."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function SIR!(du,u,p,t)\n"," S,I,R = u\n"," β,γ = p\n"," du[1] = dS = -β*S*I\n"," du[2] = dI = +β*S*I - γ*I\n"," du[3] = dR = +γ*I\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations\n","u₀ = [0.99; 0.01; 0]\n","tspan = (0.0,15.0)\n","p = (2.0,0.4)\n","problem = ODEProblem(SIR!,u₀,tspan,p)\n","solution = solve(problem)\n","plot(solution,labels=[\"susceptible\" \"infected\" \"recovered\"])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**12.14. Duffing equation.** We'll use a high-order, explicit ODE solver for the Duffing equation."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations, Plots\n","function duffing!(dx,x,γ,t) \n"," dx[1] = x[2]\n"," dx[2] = -γ*x[2]+x[1]-x[1]^3+0.3*cos(t)\n","end\n","problem = ODEProblem(duffing!,[1.0,0.0],(0.0,200.0),0.37)\n","solution = solve(problem,Vern7())\n","plot(solution, idxs=(1,2))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function solveIVP(f,u0,tspan)\n"," sol = solve(ODEProblem(f,u0,tspan))\n"," return(sol.u[end][1]) \n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using DifferentialEquations, Roots\n","airy(y,p,x) = [y[2];x*y[1]]\n","domain = (-12.0,0.0); bc = (1.0,1.0); guess = 5\n","shoot_airy = (guess -> solveIVP(airy,[bc[1];guess],domain)-bc[2])\n","v = find_zero(shoot_airy, guess)"]},{"cell_type": "markdown", "metadata": {}, "source": ["Once we have our second initial value, we can plot the solution:"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["sol = solve(ODEProblem(airy,[bc[1],v],domain))\n","plot(sol)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**13.4. An unconditionally unstable method.** Let's generate the coefficients for multistep scheme given by the stencil:
\"unstable
We'll use the function [`multistepcoefficients`](#multistepcoefficients) introduced earlier. Then we'll plot $r(\\lambda k)$ along the real axis. The method is unstable wherever $|r|>1$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m, n = [0 1 2 3 4], [1]\n","a, b = zeros(maximum(m)+1), zeros(maximum(n)+1) \n","a[m.+1], b[n.+1] = multistepcoefficients(m,n)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["λk = r -> (a ⋅ r.^-(1:length(a))) ./ (b ⋅ r.^-(1:length(b)))\n","r = LinRange(0.2,6,100)\n","plot([λk.(r) λk.(-r)],[r r],xlim=[-2,2])"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["Δx = 0.01; Δt = 0.01\n","ℓ = 1; x = -ℓ:Δx:ℓ; m = length(x)\n","Uⁿ = exp.(-8*x.^2); Uⁿ⁻¹ = Uⁿ \n","ν = Δt/Δx^2; α = 0.5 + ν; γ = 0.5 - ν\n","B = ν*Tridiagonal(ones(m-1),zeros(m),ones(m-1))\n","B[1,2] *=2; B[end,end-1] *=2\n","@gif for i = 1:300\n"," global Uⁿ, Uⁿ⁻¹ = (B*Uⁿ+γ*Uⁿ⁻¹)/α, Uⁿ\n"," plot(x,Uⁿ,ylim=(0,1),label=:none,fill=(0, 0.3, :red))\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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)$."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function schroedinger(m,n,ε)\n"," x = LinRange(-4,4,m); Δx = x[2]-x[1]; Δt = 2π/n; V = x.^2/2\n"," ψ = exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)\n"," diags = 0.5im*ε*[1 -2 1]/Δx^2 .- im/ε*[0 1 0].*V\n"," D = Tridiagonal(diags[2:end,1], diags[:,2], diags[1:end-1,3])\n"," D[1,2] *= 2; D[end,end-1] *= 2\n"," A = I + 0.5Δt*D \n"," B = I - 0.5Δt*D \n"," for i ∈ 1:n \n"," ψ = B\\(A*ψ)\n"," end\n"," return ψ\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ε = 0.3; m = 20000; eₓ=[]; eₜ=[]\n","N = floor.(Int,exp10.(LinRange(2,3.7,6)))\n","x = LinRange(-4,4,m)\n","ψₘ = -exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)\n","for n ∈ N\n"," x = LinRange(-4,4,n)\n"," ψₙ = -exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)\n"," append!(eₜ,norm(ψₘ - schroedinger(m,n,ε))/m)\n"," append!(eₓ,norm(ψₙ - schroedinger(n,m,ε))/n)\n","end\n","plot(2π./N,eₜ,shape=:circle, xaxis=:log, yaxis=:log)\n","plot!(8 ./N,eₓ,shape=:circle)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["t = 0.5; n=100; m = 100\n","r = LinRange(0,2,m); Δr = r[2]-r[1]; Δt = t/n;\n","u₀ = @. tanh(32(1 - r)); u = u₀\n","d = @. [1 -2 1]/Δr^2 + (1/r)*[-1 0 1]/2Δr\n","D = Tridiagonal(d[2:end,1],d[:,2],d[1:end-1,3])\n","D[1,1:2] = [-4 4]/Δr^2; D[end,end-1:end] = [2 -2]/Δr^2 \n","A = I - 0.5Δt*D \n","B = I + 0.5Δt*D \n","for i = 1:n\n"," u = A\\(B*u)\n","end\n","plot(r, u,label=:none,fill=(-1, 0.3, :red))"]},{"cell_type": "markdown", "metadata": {}, "source": ["The following is a slower but high-order alternative to the Crank–Nicolson method:"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Sundials\n","problem = ODEProblem((u,p,t)->D*u,u₀,(0,t))\n","method = CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)\n","solution = solve(problem,method)\n","plot(r, solution(t),label=:none,fill=(-1, 0.3, :red))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["logitspace(x,n,p) = x*atanh.(LinRange(-p,p,n))/atanh(p)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function laplacian(x)\n"," Δx = diff(x); Δx₁ = Δx[1:end-1]; Δx₂ = Δx[2:end]\n"," d₋ = @. 2/[Δx₁*(Δx₁+Δx₂); Δx₁[1]^2]\n"," d₀ = @. -2/[Δx₂[end].^2; Δx₁*Δx₂; Δx₁[1].^2] \n"," d₊ = @. 2/[Δx₂[end].^2; Δx₂*(Δx₁+Δx₂) ]\n"," Tridiagonal(d₋,d₀,d₊)\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["function heat_equation(x,t,n,u)\n"," Δt = t/n\n"," D² = laplacian(x)\n"," A = I - 0.5Δt*D²\n"," B = I + 0.5Δt*D²\n"," for i ∈ 1:n\n"," u = A\\(B*u)\n"," end\n"," return u\n","end"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["ϕ = (x,t,s) -> exp.(-s*x.^2/(1+4*s*t))/sqrt(1+4*s*t)\n","m = 40; t = 15\n","x = LinRange(-20,20,m)\n","plot(x,ϕ(x,t,10),label=\"exact\",width=3)\n","u₁ = heat_equation(x,t,m,ϕ(x,0,10))\n","plot!(x,u₁,shape=:circle,label=\"equal\")\n","x = logitspace(20,m,0.999)\n","u₂ = heat_equation(x,t,m,ϕ(x,0,10))\n","plot!(x,u₂,shape=:circle,label=\"logit\")"]},{"cell_type": "markdown", "metadata": {}, "source": ["**13.10. Allen–Cahn equation.** We'll solve the Allen–Cahn equation using Strang splitting. "]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["L = 16; m = 400; Δx = L/m\n","T = 4; n = 1600; Δt = T/n\n","x = LinRange(-L/2,L/2,m)'\n","u = @. tanh(x^4 - 16*(2*x^2-x'^2)) \n","D = Tridiagonal(ones(m-1),-2ones(m),ones(m-1))/Δx^2\n","D[1,2] *= 2; D[end,end-1] *= 2\n","A = I + 0.5Δt*D\n","B = I - 0.5Δt*D \n","f = (u,Δt) -> @. u/sqrt(u^2 - (u^2-1)*exp(-50*Δt))\n","u = f(u,Δt/2)\n","anim = Animation()\n","for i = 1:n\n"," (i%10)==1 && (plot(Gray.(u),border=:none); frame(anim))\n"," u = (B\\(A*(B\\(A*u))'))'\n"," (i @. t<2 ? \n"," (0 u.^2/2; df = u -> u\n","u = (x.>=0).*(x.<=1) \n","anim = @animate for i = 1:n \n"," α = 0.5*max.(abs.(df(u[j])),abs.(df(u[j.+1])))\n"," F = (f(u[j])+f(u[j.+1]))/2 - α.*(u[j.+1]-u[j])\n"," global u -= Δt/Δx*[0;diff(F);0]\n"," plot(x,u, fill = (0, 0.3, :blue))\n"," plot!(x,U(x,i*Δt), fill = (0, 0.3, :red)) \n"," plot!(legend=:none, ylim=[0,1])\n","end\n","gif(anim, \"burgers.gif\", fps = 15)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**14.8. Dam break problem.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["δ = u -> diff(u,dims=1)\n","ϕ = t -> (abs(t)+t)./(1+abs(t))\n","fixnan(u) = isnan(u)||isinf(u) ? 0 : u\n","θ = δu -> fixnan.(δu[1:end-1,:]./δu[2:end,:])\n","∂ₓ(u) = (δu=δ(u);[[0 0];δu[2:end,:].*ϕ.(θ(δu));[0 0]])\n","F = u -> [u[:,1].*u[:,2] u[:,1].*u[:,2].^2+0.5u[:,1].^2]"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m = 100; x = LinRange(-.5,.5,m); Δx = x[2]-x[1]\n","T = 0.4; n = ceil(T/(Δx/2)); Δt = (T/n)/2; \n","U = [0.8*(x.<0).+0.2 zero(x)] \n","Uⱼ = view(U,1:m-1,:); Uⱼ₊₁ = view(U,2:m,:)\n","anim = @animate for i = 1:n\n"," U° = U-0.5*Δt/Δx*∂ₓ(F(U))\n"," Uⱼ₊₁ .= (Uⱼ+Uⱼ₊₁)/2 - δ(∂ₓ(U))/8 - Δt/Δx*δ(F(U°)) \n"," U° = U-0.5*Δt/Δx*∂ₓ(F(U))\n"," Uⱼ .= (Uⱼ+Uⱼ₊₁)/2 - δ(∂ₓ(U))/8 - Δt/Δx*δ(F(U°))\n"," plot(x,U[:,1], fill = (0, 0.3, :blue))\n"," plot!(legend=:none, ylim=[0,1])\n","end\n","gif(anim, \"dambreak.gif\", fps = 15)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**15.1. Finite element method.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m=10; x=LinRange(0,1,m); h=x[2]-x[1]\n","α = fill(2/h-2h/3,m); α[[1,m]] /= 2; β = fill(-1/h-h/6,m-1)\n","A = SymTridiagonal(α,β)\n","b = [-2h^3/3; -4h^3/3 .- 8h*x[2:m-1].^2; -4h+8h^2/3-2h^3/3+1]\n","u = A\\b\n","s = -16 .+ 8x.^2 .+ 15csc(1).*cos.(x)\n","plot(x,s,marker=:o,alpha=0.5); plot!(x,u,marker=:o,alpha=0.5)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**15.2. Finite element method.**"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["m = 8; x = LinRange(0,1,m+2); h = x[2]-x[1]\n","σ = (a,b,c) -> Tridiagonal(fill(a,m-1),fill(b,m),fill(c,m-1))/h^3\n","M = [σ(-12,24,-12) σ(-6,0,6);σ(6,0,-6) σ(2,8,2)];\n","b = [384h*ones(m);zeros(m)]\n","u = M\\b\n","s = 16*(x.^4 - 2x.^3 + x.^2)\n","plot(x,[s [0;u[1:m];0]],marker=:o,alpha=0.5)\n","s = 16*(x.^4 - 2x.^3 + x.^2)\n","plot(x,[s [0;u[1:m];0]],marker=:o,alpha=0.5)"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, DifferentialEquations\n","m = 128; x = (1:m)/m*(2π) .- π\n","ξ = im*[0:(m÷2); (-m÷2+1):-1]\n","f = (u,p,t) -> -real(ifft(ξ.*fft(0.5u.^2)))\n","u = solve(ODEProblem(f,exp.(-x.^2),(0.0,2.0)),DP5());"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["@gif for t=0:.02:2\n"," plot(x,u(t),ylim=(-0.2,1.2))\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**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."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, DifferentialEquations\n","ϕ = (x,x₀,c) -> 0.5c*sech(sqrt(c)/2*(x-x₀))^2\n","L = 30; T = 2.0; m = 256\n","x = (1:m)/m*L .- L/2\n","u₀ = ϕ.(x,-4,4) + ϕ.(x,-9,9)\n","iξ = im*[0:(m÷2);(-m÷2+1):-1]*2π/L\n","F = plan_fft(u₀)"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["G = t -> exp.(-iξ.^3*t)\n","f = (w,_,t) -> -G(t) .\\ (3iξ .* (F * (F \\ (G(t).*w) ).^2) )"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["w = solve(ODEProblem(f,F*u₀,(0,T)),DP5())\n","u = t -> real(F\\(w(t).*G(t)))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["@gif for t=0:.01:2\n"," plot(x,u(t),ylim=(0,5),legend=:none)\n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**16.4. Swift–Hohenberg equation.** We'll use Strang splitting to solve the Swift–Hohenberg equation."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, Images\n","ϵ = 1; m = 256; ℓ = 100; n = 2000; Δt=100/n\n","U = (rand(m,m).>.5) .- 0.5\n","ξ = [0:(m÷2);(-m÷2+1):-1]*2π/ℓ\n","D²= -ξ.^2 .- (ξ.^2)'\n","E = exp.(-(D².+1).^2*Δt)\n","F = plan_fft(U)\n","f = U -> U./sqrt.(U.^2/ϵ + exp(-Δt*ϵ)*(1 .- U.^2/ϵ))\n","for i=1:600\n"," U = f(F\\(E.*(F*f(U))))\n","end\n","Gray.(clamp01.((real(U).+1)/2))"]},{"cell_type": "markdown", "metadata": {}, "source": ["**16.5. Cahn–Hilliard equation.** We'll solve the two-dimensional Cahn–Hilliard equation using a fourth-order ETDRK scheme."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, DifferentialEquations, SciMLOperators\n","ℓ = 8; T = 8.0; m = 256; ε² = 0.01; α = 1\n","ξ = [0:(m÷2);(-m÷2+1):-1]*2π/ℓ\n","Δ = -(ξ.^2 .+ ξ'.^2); Δ² = Δ.^2\n","u₀ = randn(m,m)\n","F = plan_fft(u₀)\n","L = DiagonalOperator((-ε²*Δ²+α*Δ)[:])\n","N = (u,_,_) -> (v = F\\reshape(u,m,m); Δ.*(F*(v.^3-(1+α)*v)))[:]\n","problem = SplitODEProblem(L,N,(F*u₀)[:],(0,T))\n","solution = solve(problem,ETDRK4(),dt = 0.1)\n","u = t -> real(F\\reshape(solution(t),m,m))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Interact, ColorSchemes, Images\n","@manipulate for t in slider(0:0.1:8; value=0, label=\"time\")\n"," get(colorschemes[:binary], u(t), :extrema) \n","end"]},{"cell_type": "markdown", "metadata": {}, "source": ["**16.6. Kuramoto–Sivashinsky equation.** We'll use a fourth-order ETDRK method to solve the two-dimensional KSE."]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using FFTW, DifferentialEquations, SciMLOperators\n","ℓ = 50; T = 200.0; n = 128; u₀ = randn(n,n)\n","x = LinRange(-ℓ/2,ℓ/2,n+1)[2:end]\n","ξ = [0:(n÷2);(-n÷2+1):-1]*2π/ℓ\n","Δ = -(ξ.^2 .+ ξ'.^2); Δ² = Δ.^2\n","F = plan_fft(u₀)\n","L = DiagonalOperator((-Δ-Δ²)[:])\n","N = (u,_,_) -> \n"," ( v = reshape(u,n,n); \n"," v = -0.5*abs.((F\\(im*ξ.*v)).^2 + (F\\(im*ξ'.*v)).^2);\n"," w = (F*v)[:]; w[1] = 0.0; return w )\n","problem = SplitODEProblem(L,N,(F*u₀)[:],(0,T))\n","solution = solve(problem, ETDRK4(), dt=0.05, saveat=0.5);\n","u = t -> real(F\\reshape(solution(t),n,n))"]}, {"cell_type": "code", "execution_count": null,"metadata": {},"outputs": [],"source": ["using Interact, ColorSchemes, Images\n","@manipulate for t in slider(0:0.5:T; value=0, label=\"time\")\n"," imresize(get(colorschemes[:magma], -u(t), :extrema),(256,256)) \n","end"]}], "metadata": { "kernelspec": { "display_name": "Julia 1.9.2", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.2" } }, "nbformat": 4, "nbformat_minor": 2 }