pluto_nb_SIAMFANL.jl

### 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