Skip to content

Instantly share code, notes, and snippets.

@briochemc
Last active April 30, 2021 03:17
Show Gist options
  • Save briochemc/7bb763bb64a9639d34ad7c5136356106 to your computer and use it in GitHub Desktop.
Save briochemc/7bb763bb64a9639d34ad7c5136356106 to your computer and use it in GitHub Desktop.
### A Pluto.jl notebook ###
# v0.14.4
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
el
end
end
# ╔═╡ fa21dace-a886-11eb-0a3e-8dacc8154566
begin
# Setting up the environment and packages
import Pkg
Pkg.activate(mktempdir())
Pkg.add("SIAMFANLEquations")
# This is a fork where I added a Project file
Pkg.add(url="https://github.com/briochemc/NotebookSIAMFANL.git")
Pkg.add("BandedMatrices")
Pkg.add("BenchmarkTools")
Pkg.add("AbstractFFTs")
Pkg.add("LaTeXStrings")
Pkg.add("FFTW")
Pkg.add("Plots")
Pkg.add("DataFrames")
Pkg.add("PlutoUI")
# using stdlibs
using LinearAlgebra
using SuiteSparse
using SparseArrays
# using packages
using SIAMFANLEquations
using SIAMFANLEquations.TestProblems
using SIAMFANLEquations.Examples
using NotebookSIAMFANL
using BandedMatrices
using BenchmarkTools
using AbstractFFTs
using LaTeXStrings
using FFTW
using Plots
using DataFrames
using PlutoUI
end
# ╔═╡ 5765b4e0-060a-4d77-a9a6-7842e11cf638
md"# A Pluto notebook using the [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl) solvers"
# ╔═╡ e2bb2660-63f0-45d3-b650-070eb1179a61
md"## Chapter 1"
# ╔═╡ 14b6630f-306e-4a8e-bd3f-36f1537a62e6
md"### Solver 1, `nsolc`"
# ╔═╡ 276053c8-1321-48af-ad68-ac4e47d0a5e8
md"#### Example 1, $f(x) = \tan(x) - x$"
# ╔═╡ 5dd40812-3f61-4d46-8947-6edb665dfd34
f_Ch1_alg1_ex1(x) = tan(x) - x
# ╔═╡ 4d86b416-bd88-45f8-9b90-5612013e2919
md"x₀ = $(@bind x0_Ch1_alg1_ex1 Slider(4.0:0.0001:4.7, show_value=true, default=4.5))"
# ╔═╡ e8a206c4-82f2-4a5d-9101-c3ad22dded7f
sol_Ch1_alg1_ex1 = nsolsc(f_Ch1_alg1_ex1, x0_Ch1_alg1_ex1; maxit=14, rtol=1e-17, atol=1e-17, printerr=false, keepsolhist=true)
# ╔═╡ 2d7bbbc0-371d-48c3-aeb2-d3b8bab163eb
md"#### Example 2, $f(x) = \arctan(x)$"
# ╔═╡ 50c8db3f-0820-424a-9519-e21f7f9ee7cb
f_Ch1_alg1_ex2 = atan
# ╔═╡ fd5dabb8-eb6e-4bfd-ac17-c67687826061
md"x₀ = $(@bind x0_Ch1_alg1_ex2 Slider(range(-5, 5, length=1000), show_value=true, default=1.0))"
# ╔═╡ e52d62b9-62ed-40db-aae6-c668b0862d5e
sol_Ch1_alg1_ex2 = nsolsc(f_Ch1_alg1_ex2, x0_Ch1_alg1_ex2; maxit=5, atol=1.e-12, rtol=1.e-12, keepsolhist=true)
# ╔═╡ 294b2220-abf7-41fb-a8c1-cd8780b3ccfb
md"#### Example 3, $f(x) = x^2 - 4$ and $f'(x) = 2x$"
# ╔═╡ d4c0c5d3-d8d6-4af7-9cb2-236e4e00e961
begin
f_Ch1_alg1_ex3(x) = x^2 - 4
f_Ch1_alg1_ex3′(x) = 2x
end
# ╔═╡ 171ead0d-9f69-4cd5-991a-25662a9f49f9
md"x₀ = $(@bind x0_Ch1_alg1_ex3 Slider(range(-5, 5, length=1000), show_value=true, default=1.0))"
# ╔═╡ 564cbddd-b8f3-4f3b-af26-771753f6d575
sol_Ch1_alg1_ex3 = nsolsc(f_Ch1_alg1_ex3, x0_Ch1_alg1_ex3, f_Ch1_alg1_ex3′; maxit=5, atol=1.e-9, rtol=1.e-9, keepsolhist=true)
# ╔═╡ 907202a7-2850-40b7-86b0-b023abb4dd4a
md"### Solver 2, `secant`"
# ╔═╡ ca8564f7-9b7f-400e-9da6-32f7fe5f2f30
f_Ch1_alg2_ex1 = atan
# ╔═╡ 3d983103-8f87-4d25-baec-13488fe21d05
md"x₀ = $(@bind x0_Ch1_alg2_ex1 Slider(range(-5, 5, length=1000), show_value=true, default=1.0))"
# ╔═╡ 6a9c94ec-83b8-4f24-8e38-158d6ab271a8
sol_Ch1_alg2_ex1 = secant(f_Ch1_alg2_ex1, x0_Ch1_alg2_ex1; maxit=6, atol=1.e-12, rtol=1.e-12, keepsolhist=true)
# ╔═╡ 9e2d52b8-2c93-4472-b4a3-69ec0d8cf958
md"### Comparing solvers"
# ╔═╡ 5276f299-b45a-4322-a413-e857dee5ece6
md"x₀ = $(@bind x0_Ch1_alg3_ex1 Slider(range(-5, 5, length=10001), show_value=true, default=1.0))"
# ╔═╡ 3788fcfc-ed24-48c7-8909-1d69ff27f9ee
begin
f_Ch1_alg3_ex1 = atan
opts = (rtol=1.e-12, atol=1.e-12)
sol_Ch1_alg3_ex1 = Dict(
:newton => nsolsc(f_Ch1_alg3_ex1, x0_Ch1_alg3_ex1; opts...),
:secant => secant(f_Ch1_alg3_ex1, x0_Ch1_alg3_ex1; opts...),
:chord => nsolsc(f_Ch1_alg3_ex1, x0_Ch1_alg3_ex1; opts..., solver="chord")
)
end
# ╔═╡ c55fa796-0173-4be1-afc1-6c86e2a659a0
md"## Chapter 2"
# ╔═╡ 80ef44e2-4176-4d2a-82ab-821b67826ad5
md"## Extra functions"
# ╔═╡ 6f2b0404-939c-4699-bac8-7fd6856939c0
# Convenience function to extend the xlims and ylims window when plotting
function extended_window(lims)
min, max = lims
Δ = max - min
min - Δ/10, max + Δ/10
end
# ╔═╡ 546d2e8e-c76b-489a-9e19-6e7f24133648
# Convenience function to extend the xlims and ylims window when plotting
function extended_window(lims_list::Vector)
extrema(collect(extended_window(lims)) for lims in lims_list)
end
# ╔═╡ 45eec6b7-c7ca-4ea9-b40e-9d28240caf1d
function myplot(f, sol; kwargs...)
xlims = extended_window(extrema(sol.solhist))
ylims = extended_window(extrema(f.(sol.solhist)))
xs = range(get(kwargs, :xlims, xlims)..., length=1000) # if given use xlims kwarg
p = plot(xs, f, lab=L"f(x)"; ylims, xlims)
hline!(p, [0], lab="", ls=:dash, c=:black)
plot!(p, sol.solhist', f.(sol.solhist'), st=:steppost, leg=:outerright, lab="$(length(sol.solhist)) iterations", xlab=L"x", ylab=L"f(x)", m=:o, ms=2, c=2; kwargs...)
p
end
# ╔═╡ 8c61b8b7-2235-480d-b83e-eac678bc4fa6
function smiley(sol)
# errcode = 0 if the iteration succeeded
# = -1 if the initial iterate satisifies the termination criteria
# = 10 if no convergence after maxit iterations
# = 1 if the line search failed
e = sol.errcode
return (e == 0) ? "✔" : ((e == 1) ? "✖ (line search)" : ((e== 10) ? "✖ (maxit)" : "✔ (x₀ satisfies)"))
end
# ╔═╡ aba5b16f-f304-4039-80b1-28013164bb77
function myplot(f, sols::Dict; kwargs...)
xlims = extended_window([extrema(sol.solhist) for (k,sol) in sols])
ylims = extended_window([extrema(f.(sol.solhist)) for (k,sol) in sols])
xs = range(get(kwargs, :xlims, xlims)..., length=1000) # if given use xlims kwarg
p = plot(xs, f, lab=L"f(x)"; ylims, xlims)
hline!(p, [0], lab="", ls=:dash, c=:black)
c = 2
for (k,sol) in sols
plot!(p, sol.solhist', f.(sol.solhist'), ls=c, c=c, st=:steppost, leg=:outerbottom, lab="$k: $(smiley(sol)) $(length(sol.solhist)) iterations", xlab=L"x", ylab=L"f(x)", α=0.9, m=:o, ms=2; kwargs...)
c += 1
end
p
end
# ╔═╡ 33da3a19-c5bc-46a0-9b70-6ace42d63b2e
myplot(f_Ch1_alg1_ex1, sol_Ch1_alg1_ex1; title=L"\textrm{Solving }f(x) = \tan(x)-x=0")
# ╔═╡ a00014f9-f63e-40b7-958d-a8e1d038b840
myplot(f_Ch1_alg1_ex2, sol_Ch1_alg1_ex2, title=L"\textrm{Solving }f(x) = \arctan(x) = 0")
# ╔═╡ 130f2e88-f6aa-4320-a2f9-a131760c2e1e
myplot(f_Ch1_alg1_ex3, sol_Ch1_alg1_ex3; xlims=(-5,5), ylims=(-5,5), title=L"\textrm{Solving }f(x)=x^2-4 = 0 \textrm{ and providing }f'(x) = 2x")
# ╔═╡ 5637659b-3b03-4c66-927e-44bc59ee9a19
myplot(f_Ch1_alg2_ex1, sol_Ch1_alg2_ex1, title=L"\textrm{Sovling } f(x) = tan(x) \textrm{ with the secant solver}", xlims=(-5,5), ylims=(-π/2,π/2))
# ╔═╡ b7a8d69a-3fb8-49e6-abf8-7badc9249a2d
myplot(f_Ch1_alg3_ex1, sol_Ch1_alg3_ex1, title="Comparing solvers", xlims=(-5,5), ylims=(-π/2, π/2))
# ╔═╡ Cell order:
# ╟─fa21dace-a886-11eb-0a3e-8dacc8154566
# ╟─5765b4e0-060a-4d77-a9a6-7842e11cf638
# ╟─e2bb2660-63f0-45d3-b650-070eb1179a61
# ╟─14b6630f-306e-4a8e-bd3f-36f1537a62e6
# ╟─276053c8-1321-48af-ad68-ac4e47d0a5e8
# ╠═5dd40812-3f61-4d46-8947-6edb665dfd34
# ╠═e8a206c4-82f2-4a5d-9101-c3ad22dded7f
# ╠═4d86b416-bd88-45f8-9b90-5612013e2919
# ╠═33da3a19-c5bc-46a0-9b70-6ace42d63b2e
# ╟─2d7bbbc0-371d-48c3-aeb2-d3b8bab163eb
# ╠═50c8db3f-0820-424a-9519-e21f7f9ee7cb
# ╠═e52d62b9-62ed-40db-aae6-c668b0862d5e
# ╠═fd5dabb8-eb6e-4bfd-ac17-c67687826061
# ╟─a00014f9-f63e-40b7-958d-a8e1d038b840
# ╟─294b2220-abf7-41fb-a8c1-cd8780b3ccfb
# ╠═d4c0c5d3-d8d6-4af7-9cb2-236e4e00e961
# ╠═564cbddd-b8f3-4f3b-af26-771753f6d575
# ╟─171ead0d-9f69-4cd5-991a-25662a9f49f9
# ╟─130f2e88-f6aa-4320-a2f9-a131760c2e1e
# ╟─907202a7-2850-40b7-86b0-b023abb4dd4a
# ╠═ca8564f7-9b7f-400e-9da6-32f7fe5f2f30
# ╠═6a9c94ec-83b8-4f24-8e38-158d6ab271a8
# ╟─3d983103-8f87-4d25-baec-13488fe21d05
# ╠═5637659b-3b03-4c66-927e-44bc59ee9a19
# ╟─9e2d52b8-2c93-4472-b4a3-69ec0d8cf958
# ╠═3788fcfc-ed24-48c7-8909-1d69ff27f9ee
# ╠═5276f299-b45a-4322-a413-e857dee5ece6
# ╠═b7a8d69a-3fb8-49e6-abf8-7badc9249a2d
# ╟─c55fa796-0173-4be1-afc1-6c86e2a659a0
# ╟─80ef44e2-4176-4d2a-82ab-821b67826ad5
# ╠═6f2b0404-939c-4699-bac8-7fd6856939c0
# ╠═546d2e8e-c76b-489a-9e19-6e7f24133648
# ╠═45eec6b7-c7ca-4ea9-b40e-9d28240caf1d
# ╠═aba5b16f-f304-4039-80b1-28013164bb77
# ╠═8c61b8b7-2235-480d-b83e-eac678bc4fa6
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment