daviddoji/CalculusWithJuliaNotes.jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

orthogonal; work around plotly()

Browse Source
This commit is contained in:
jverzani
2025-06-27 14:29:32 -04:00
parent c8f9fd4995
commit 580e87ccb2
13 changed files with 853 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
quarto/_quarto.yml
View File

@@ -83,6 +83,7 @@ book:
- integrals/volumes_slice.qmd
- integrals/arc_length.qmd
- integrals/surface_area.qmd
- integrals/orthogonal_polynomials.qmd
- integrals/twelve-qs.qmd
- part: ODEs.qmd
@@ -115,7 +116,7 @@ book:
chapters:
- alternatives/symbolics.qmd
- alternatives/SciML.qmd
# - alternatives/interval_arithmetic.qmd
#- alternatives/interval_arithmetic.qmd
- alternatives/plotly_plotting.qmd
- alternatives/makie_plotting.qmd
@@ -157,5 +158,7 @@ format:
execute:
error: false
# freeze: false
freeze: auto
freeze: false
# freeze: auto
# cache: false
# enabled: true

38
quarto/adjust_plotly.jl
View File

@@ -5,19 +5,39 @@
# This little script just adds a line *before* the require call
# which seems to make it all work. The line number 83 might change.
f = "_book/alternatives/plotly_plotting.html"
lineno = 88
#alternatives/plotly_plotting.html
function _add_plotly(f)
lineno = 117
str = """
<script src="https://cdn.plot.ly/plotly-2.11.0.min.js"></script>
"""
r = readlines(f)
open(f, "w") do io
for (i,l) ∈ enumerate(r)
i == lineno && println(io, str)
println(io, l)
r = readlines(f)
open(f, "w") do io
for (i,l) ∈ enumerate(r)
i == lineno && println(io, str)
println(io, l)
end
end
end
function (@main)(args...)
for (root, dirs, files) in walkdir("_book")
for fᵢ ∈ files
f = joinpath(root, fᵢ)
if endswith(f, ".html")
_add_plotly(f)
end
end
end
#f = "_book/integrals/center_of_mass.html"
#_add_plotly(f)
return 1
end
["ODEs", "alternatives", "derivatives", "differentiable_vector_calculus", "integral_vector_calculus", "integrals", "limits", "misc", "precalc", "site_libs"]

4
quarto/alternatives.qmd
View File

@@ -6,8 +6,8 @@ These notes use a particular selection of packages. This selection could have be
* The finding of zeros of scalar-valued, univariate functions is done with `Roots`. The [NonlinearSolve](./alternatives/SciML.html#nonlinearsolve) package provides an alternative for univariate and multi-variate functions.
* The finding of minima and maxima was done mirroring the framework of a typical calculus class; the [Optimization](./alternatives/SciML.html#optimization-optimization.jl) provides an alternative.
* The finding of minima and maxima was done mirroring the framework of a typical calculus class; the [Optimization](./alternatives/SciML.html#optimization-optimization.jl) package provides an alternative.
* The computation of numeric approximations for definite integrals is computed with the `QuadGK` and `HCubature` packages. The [Integrals](./alternatives/SciML.html#integration-integrals.jl) package provides a unified interface for numeric to these two packages, among others.
* The computation of numeric approximations for definite integrals is computed with the `QuadGK` and `HCubature` packages. The [Integrals](./alternatives/SciML.html#integration-integrals.jl) package provides a unified interface for numeric integration, including these two packages, among others.
* Plotting was done using the popular `Plots` package. The [Makie](./alternatives/makie_plotting.html) package provides a very powerful alternative. Whereas the [PlotlyLight](./alternatives/plotly_plotting.html) package provides a light-weight alternative using an open-source JavaScript library.

7
quarto/alternatives/Project.toml
View File

@@ -5,11 +5,13 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
Implicit3DPlotting = "d997a800-832a-4a4c-b340-7dddf3c1ad50"
Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
Meshing = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
@@ -19,6 +21,7 @@ PlotlyKaleido = "f2990250-8cf9-495f-b13a-cce12b45703c"
PlotlyLight = "ca7969ec-10b3-423e-8d99-40f33abb42bf"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
QuizQuestions = "612c44de-1021-4a21-84fb-7261cf5eb2d4"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -26,3 +29,5 @@ SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
SymbolicLimits = "19f23fe9-fdab-4a78-91af-e7b7767979c3"
SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
TextWrap = "b718987f-49a8-5099-9789-dcd902bef87d"

2
quarto/alternatives/SciML.qmd
View File

@@ -250,7 +250,7 @@ With the system defined, we can pass this to `NonlinearProblem`, as was done wit
```{julia}
prob = NonlinearProblem(ns, [1.0], [α => 1.0])
prob = NonlinearProblem(mtkcompile(ns), [1.0], Dict(α => 1.0))
```
The problem is solved as before:

80
quarto/alternatives/makie_plotting.qmd
View File

@@ -1,6 +1,5 @@
# Calculus plots with Makie
{{< include ../_common_code.qmd >}}
The [Makie.jl webpage](https://github.com/JuliaPlots/Makie.jl) says
@@ -36,8 +35,7 @@ using GLMakie
import LinearAlgebra: norm
```
The `Makie` developers have workarounds for the delayed time to first plot, but without utilizing these the time to load the package is lengthy.
The package load time as of recent version of `Makie` is quite reasonable for a complicated project. (The time to first plot is under 3 seconds on a typical machine.)
## Points (`scatter`)
@@ -158,7 +156,8 @@ A point is drawn with a "marker" with a certain size and color. These attributes
```{julia}
scatter(xs, ys;
marker=[:x,:cross, :circle], markersize=25,
marker=[:x,:cross, :circle],
markersize=25,
color=:blue)
```
@@ -176,7 +175,7 @@ A single value will be repeated. A vector of values of a matching size will spec
## Curves
The curves of calculus are lines. The `lines` command of `Makie` will render a curve by connecting a series of points with straight-line segments. By taking a sufficient number of points the connect-the-dot figure can appear curved.
A visualization of a curve in calculus is comprised of line segments. The `lines` command of `Makie` will render a curve by connecting a series of points with straight-line segments. By taking a sufficient number of points the connect-the-dot figure can appear curved.
### Plots of univariate functions
@@ -304,6 +303,7 @@ current_figure()
### Text (`annotations`)
XXX FIX ME XXX
Text can be placed at a point, as a marker is. To place text, the desired text and a position need to be specified along with any adjustments to the default attributes.
@@ -315,25 +315,43 @@ For example:
xs = 1:5
pts = Point2.(xs, xs)
scatter(pts)
annotations!("Point " .* string.(xs), pts;
fontsize = 50 .- 2*xs,
rotation = 2pi ./ xs)
annotation!(pts;
text = "Point " .* string.(xs),
fontsize = 30 .- 5*xs)
current_figure()
```
The graphic shows that `fontsize` adjusts the displayed size and `rotation` adjusts the orientation. (The graphic also shows a need to manually override the limits of the `y` axis, as the `Point 5` is chopped off; the `ylims!` function to do so will be shown later.)
The graphic shows that `fontsize` adjusts the displayed size.
Attributes for `text`, among many others, include:
* `align` Specify the text alignment through `(:pos, :pos)`, where `:pos` can be `:left`, `:center`, or `:right`.
* `rotation` to indicate how the text is to be rotated
* `fontsize` the font point size for the text
* `font` to indicate the desired font
Annotations with an arrow can be useful to highlight a feature of a graph. This example is modified from the documentation and utilizes some interval functions to draw an arrow with an arc:
```{julia}
g(x) = cos(6x) * exp(x)
xs = 0:0.01:4
_, ax, _ = lines(xs, g.(xs); axis = (; xgridvisible = false, ygridvisible = false))
annotation!(ax, 1, 20, 2.1, g(2.1),
text = "A relative maximum",
path = Ann.Paths.Arc(0.3),
style = Ann.Styles.LineArrow(),
labelspace = :data
)
current_figure()
```
#### Line attributes
@@ -666,6 +684,7 @@ A surface of revolution for $g(u)$ revolved about the $z$ axis can be visualized
```{julia}
g(u) = u^2 * exp(-u)
r(u,v) = (g(u)*sin(v), g(u)*cos(v), u)
us = range(0, 3, length=10)
vs = range(0, 2pi, length=10)
xs, ys, zs = parametric_grid(us, vs, r)
@@ -681,6 +700,7 @@ A torus with big radius $2$ and inner radius $1/2$ can be visualized as follows
```{julia}
r1, r2 = 2, 1/2
r(u,v) = ((r1 + r2*cos(v))*cos(u), (r1 + r2*cos(v))*sin(u), r2*sin(v))
us = vs = range(0, 2pi, length=25)
xs, ys, zs = parametric_grid(us, vs, r)
@@ -696,6 +716,7 @@ A Möbius strip can be produced with:
ws = range(-1/4, 1/4, length=8)
thetas = range(0, 2pi, length=30)
r(w, θ) = ((1+w*cos(θ/2))*cos(θ), (1+w*cos(θ/2))*sin(θ), w*sin(θ/2))
xs, ys, zs = parametric_grid(ws, thetas, r)
surface(xs, ys, zs)
@@ -865,20 +886,19 @@ end
#### Implicitly defined surfaces, $F(x,y,z)=0$
To plot the equation $F(x,y,z)=0$, for $F$ a scalar-valued function, again the implicit function theorem says that, under conditions, near any solution $(x,y,z)$, $z$ can be represented as a function of $x$ and $y$, so the graph will look likes surfaces stitched together. The `Implicit3DPlotting` package takes an approach like `ImplicitPlots` to represent these surfaces. It replaces the `Contour` package computation with a $3$-dimensional alternative provided through the `Meshing` and `GeometryBasics` packages.
To plot the equation $F(x,y,z)=0$, for $F$ a scalar-valued function, again the implicit function theorem says that, under conditions, near any solution $(x,y,z)$, $z$ can be represented as a function of $x$ and $y$, so the graph will look like surfaces stitched together.
```{julia}
using Implicit3DPlotting
```
With `Makie`, many implicitly defined surfaces can be adequately represented using `countour` with the attribute `levels=[0]`. We will illustrate this technique.
The `Implicit3DPlotting` package takes an approach like `ImplicitPlots` to represent these surfaces. It replaces the `Contour` package computation with a $3$-dimensional alternative provided through the `Meshing` and `GeometryBasics` packages. This package has a `plot_implicit_surface` function that does something similar to below. We don't illustrate it, as it *currently* doesn't work with the latest version of `Makie`.
This example, plotting an implicitly defined sphere, comes from the documentation of `Implicit3DPlotting`. The `f` to be plotted is a scalar-valued function of a vector:
To begin, we plot a sphere implicitly as a solution to $F(x,y,z) = x^2 + y^2 + z^2 - 1 = 0$>
```{julia}
f(x) = sum(x.^2) - 1
xlims = ylims = zlims = (-5, 5)
plot_implicit_surface(f; xlims, ylims, zlims)
f(x,y,z) = x^2 + y^2 + z^2 - 1
xs = ys = zs = range(-3/2, 3/2, 100)
contour(xs, ys, zs, f; levels=[0])
```
@@ -887,11 +907,13 @@ Here we visualize an intersection of a sphere with another figure:
```{julia}
r₂(x) = sum(x.^2) - 5/4 # a sphere
r₂(x) = sum(x.^2) - 2 # a sphere
r₄(x) = sum(x.^4) - 1
xlims = ylims = zlims = (-2, 2)
p = plot_implicit_surface(r₂; xlims, ylims, zlims, color=:yellow)
plot_implicit_surface!(p, r₄; xlims, ylims, zlims, color=:red)
ϕ(x,y,z) = (x,y,z)
xs = ys = zs = range(-2, 2, 100)
contour(xs, ys, zs, r₂∘ϕ; levels = [0], colormap=:RdBu)
contour!(xs, ys, zs, r₄∘ϕ; levels = [0], colormap=:viridis)
current_figure()
```
@@ -900,11 +922,12 @@ This example comes from [Wikipedia](https://en.wikipedia.org/wiki/Implicit_surfa
```{julia}
f(x,y,z) = 2y*(y^2 -3x^2)*(1-z^2) + (x^2 +y^2)^2 - (9z^2-1)*(1-z^2)
xlims = ylims = zlims = (-5/2, 5/2)
plot_implicit_surface(x -> f(x...); xlims, ylims, zlims)
xs = ys = zs = range(-5/2, 5/2, 100)
contour(xs, ys, zs, f; levels=[0], colormap=:RdBu)
```
(This figure does not render well through `contour(xs, ys, zs, f, levels=[0])`, as the hole is not shown.)
(This figure does not render well though, as the hole is not shown.)
For one last example from Wikipedia, we have the Cassini oval which "can be defined as the point set for which the *product* of the distances to $n$ given points is constant." That is:
@@ -915,8 +938,8 @@ function cassini(λ, ps = ((1,0,0), (-1, 0, 0)))
n = length(ps)
x -> prod(norm(x .- p) for p ∈ ps) - λ^n
end
xlims = ylims = zlims = (-2, 2)
plot_implicit_surface(cassini(1.05); xlims, ylims, zlims)
xs = ys = zs = range(-2, 2, 100)
contour(xs, ys, zs, cassini(0.80) ∘ ϕ; levels=[0], colormap=:RdBu)
```
## Vector fields. Visualizations of $f:R^2 \rightarrow R^2$
@@ -1064,7 +1087,7 @@ F
### Observables
The basic components of a plot in `Makie` can be updated [interactively](https://makie.juliaplots.org/stable/documentation/nodes/index.html#observables_interaction). `Makie` uses the `Observables` package which allows complicated interactions to be modeled quite naturally. In the following we give a simple example.
The basic components of a plot in `Makie` can be updated [interactively](https://makie.juliaplots.org/stable/documentation/nodes/index.html#observables_interaction). Historically `Makie` used the `Observables` package which allows complicated interactions to be modeled quite naturally. In the following we give a simple example, though newer versions of `Makie` rely on a different mechanism.
In Makie, an `Observable` is a structure that allows its value to be updated, similar to an array. When changed, observables can trigger an event. Observables can rely on other observables, so events can be cascaded.
@@ -1123,6 +1146,7 @@ end
lines!(ax, xs, f)
lines!(ax, points)
scatter!(ax, points; markersize=10)
current_figure()
```

58
quarto/alternatives/plotly_plotting.qmd
View File

@@ -73,7 +73,6 @@ The `Config` constructor (from the `EasyConfig` package loaded with `PlotlyLight
```{julia}
#| hold: true
cfg = Config()
cfg.key1.key2.key3 = "value"
cfg
@@ -89,7 +88,6 @@ A basic scatter plot of points $(x,y)$ is created as follows:
```{julia}
#| hold: true
xs = 1:5
ys = rand(5)
data = Config(x = xs,
@@ -113,7 +111,6 @@ A line plot is very similar, save for a different `mode` specification:
```{julia}
#| hold: true
xs = 1:5
ys = rand(5)
data = Config(x = xs,
@@ -134,7 +131,6 @@ The line graph plays connect-the-dots with the points specified by paired `x` an
```{julia}
#| hold: true
data = Config(
x=[0,1,nothing,3,4,5],
y = [0,1,2,3,4,5],
@@ -149,7 +145,6 @@ More than one graph or layer can appear on a plot. The `data` argument can be a
```{julia}
#| hold: true
data = [Config(x = 1:5,
y = rand(5),
type = "scatter",
@@ -177,7 +172,6 @@ For example, here we plot the graphs of both the $\sin(x)$ and $\cos(x)$ over $[
```{julia}
#| hold: true
a, b = 0, 2pi
xs, ys = PlotUtils.adapted_grid(sin, (a,b))
@@ -193,7 +187,6 @@ The values for `a` and `b` are used to generate the $x$- and $y$-values. These c
```{julia}
#| hold: true
xs, ys = PlotUtils.adapted_grid(x -> x^5 - x - 1, (0, 2)) # answer is (0,2)
p = Plot([Config(x=xs, y=ys, name="Polynomial"),
Config(x=xs, y=0 .* ys, name="x-axis", mode="lines", line=Config(width=5))]
@@ -232,7 +225,6 @@ A marker's attributes can be adjusted by values passed to the `marker` key. Labe
```{julia}
#| hold: true
data = Config(x = 1:5,
y = rand(5),
mode="markers+text",
@@ -251,40 +243,7 @@ The `text` mode specification is necessary to have text be displayed on the char
#### RGB Colors
The `ColorTypes` package is the standard `Julia` package providing an `RGB` type (among others) for specifying red-green-blue colors. To make this work with `Config` and `JSON3` requires some type-piracy (modifying `Base.string` for the `RGB` type) to get, say, `RGB(0.5, 0.5, 0.5)` to output as `"rgb(0.5, 0.5, 0.5)"`. (RGB values in JavaScript are integers between $0$ and $255$ or floating point values between $0$ and $1$.) A string with this content can be specified. Otherwise, something like the following can be used to avoid the type piracy:
```{julia}
struct rgb
r
g
b
end
PlotlyLight.JSON3.StructTypes.StructType(::Type{rgb}) = PlotlyLight.JSON3.StructTypes.StringType()
Base.string(x::rgb) = "rgb($(x.r), $(x.g), $(x.b))"
```
With these defined, red-green-blue values can be used for colors. For example to give a range of colors, we might have:
```{julia}
#| hold: true
cols = [rgb(i,i,i) for i in range(10, 245, length=5)]
sizes = [12, 16, 20, 24, 28]
data = Config(x = 1:5,
y = rand(5),
mode="markers+text",
type="scatter",
name="scatter plot",
text = ["marker $i" for i in 1:5],
textposition = "top center",
marker = Config(size=sizes, color=cols)
)
Plot(data)
```
The `opacity` key can be used to control the transparency, with a value between $0$ and $1$.
The `ColorTypes` package is the standard `Julia` package providing an `RGB` type (among others) for specifying red-green-blue colors. To make this work with `Config` and `JSON3` requires some type-piracy (modifying `Base.string` for the `RGB` type) to get, say, `RGB(0.5, 0.5, 0.5)` to output as `"rgb(0.5, 0.5, 0.5)"`. (RGB values in JavaScript are integers between $0$ and $255$ or floating point values between $0$ and $1$.) A string with this content can be specified.
#### Marker symbols
@@ -293,7 +252,6 @@ The `marker_symbol` key can be used to set a marker shape, with the basic values
```{julia}
#| hold: true
markers = ["circle", "square", "diamond", "cross", "x", "triangle", "pentagon",
"hexagram", "star", "diamond", "hourglass", "bowtie", "asterisk",
"hash", "y", "line"]
@@ -327,7 +285,6 @@ The `shape` attribute determine how the points are connected. The default is `li
```{julia}
#| hold: true
shapes = ["linear", "hv", "vh", "hvh", "vhv", "spline"]
data = [Config(x = 1:5, y = 5*(i-1) .+ [1,3,2,3,1], mode="lines+markers", type="scatter",
name=shape,
@@ -358,7 +315,6 @@ In the following, to highlight the difference between $f(x) = \cos(x)$ and $p(x)
```{julia}
#| hold: true
xs = range(-1, 1, 100)
data = [
Config(
@@ -381,7 +337,6 @@ The `toself` declaration is used below to fill in a polygon:
```{julia}
#| hold: true
data = Config(
x=[-1,1,1,-1,-1], y = [-1,1,-1,1,-1],
fill="toself",
@@ -399,7 +354,6 @@ The legend is shown when $2$ or more charts or specified, by default. This can b
```{julia}
#| hold: true
data = Config(x=1:5, y=rand(5), type="scatter", mode="markers", name="legend label")
lyt = Config(title = "Main chart title",
xaxis = Config(title="x-axis label"),
@@ -416,7 +370,6 @@ The aspect ratio of the chart can be set to be equal through the `scaleanchor` k
```{julia}
#| hold: true
ts = range(0, 2pi, length=100)
data = Config(x = sin.(ts), y = cos.(ts), mode="lines", type="scatter")
lyt = Config(title = "A circle",
@@ -434,7 +387,6 @@ Text annotations may be specified as part of the layout object. Annotations may
```{julia}
#| hold: true
data = Config(x = [0, 1], y = [0, 1], mode="markers", type="scatter")
layout = Config(title = "Annotations",
xaxis = Config(title="x",
@@ -452,7 +404,7 @@ Plot(data, layout)
The following example is more complicated use of the elements previously described. It mimics an image from [Wikipedia](https://en.wikipedia.org/wiki/List_of_trigonometric_identities) for trigonometric identities. The use of `LaTeX` does not seem to be supported through the `JavaScript` interface; unicode symbols are used instead. The `xanchor` and `yanchor` keys are used to position annotations away from the default. The `textangle` key is used to rotate text, as desired.
```{julia, hold=true}
```{julia}
alpha = pi/6
beta = pi/5
xₘ = cos(alpha)*cos(beta)
@@ -569,7 +521,6 @@ Earlier, we plotted a two dimensional circle, here we plot the related helix.
```{julia}
#| hold: true
helix(t) = [cos(t), sin(t), t]
ts = range(0, 4pi, length=200)
@@ -596,7 +547,6 @@ There is no `quiver` plot for `plotly` using JavaScript. In $2$-dimensions a tex
```{julia}
#| hold: true
helix(t) = [cos(t), sin(t), t]
helix′(t) = [-sin(t), cos(t), 1]
ts = range(0, 4pi, length=200)
@@ -642,7 +592,6 @@ A contour plot is created by the "contour" trace type. The data is prepared as a
```{julia}
#| hold: true
f(x,y) = x^2 - 2y^2
xs = range(0,2,length=25)
@@ -661,7 +610,6 @@ The same `zs` data can be achieved by broadcasting and then collecting as follow
```{julia}
#| hold: true
f(x,y) = x^2 - 2y^2
xs = range(0,2,length=25)
@@ -692,7 +640,6 @@ Surfaces defined through a scalar-valued function are drawn quite naturally, sav
```{julia}
#| hold: true
peaks(x,y) = 3 * (1-x)^2 * exp(-(x^2) - (y+1)^2) -
10*(x/5 - x^3 - y^5) * exp(-x^2-y^2) - 1/3 * exp(-(x+1)^2 - y^2)
@@ -713,7 +660,6 @@ For parametrically defined surfaces, the $x$ and $y$ values also correspond to m
```{julia}
#| hold: true
r, R = 1, 5
X(theta,phi) = [(r*cos(theta)+R)*cos(phi),
(r*cos(theta)+R)*sin(phi),

9
quarto/alternatives/symbolics.qmd
View File

@@ -601,14 +601,7 @@ det(N)
det(collect(N))
```
Similarly, with `norm`:
```{julia}
norm(v)
```
and
Similarly, with `norm`, which returns a generator unless collected:
```{julia}

16
quarto/differentiable_vector_calculus/matrix_calculus_notes.qmd
View File

@@ -355,9 +355,15 @@ The operator $f'(x)= g'(h(x)) h'(x)$ is a product of matrices.
### Computational differences with expressions from the chain rule
Of note here is the application of the chain rule to three (or more compositions):
Of note here is the application of the chain rule to three (or more compositions) where $c:R^n \rightarrow R^j$, $b:R^j \rightarrow R^k$, and $a:R^k \rightarrow R^m$:
The derivative of $f(x) = a(x) b(x) c(x)$ can be expressed as
If $f(x) = a(b(c(x)))$ then the derivative is:
$$
f'(x) = a'(b(c(x))) b'(c(x)) c'(x),
$$
which can be expressed as three matrix multiplications two ways:
$$
f' = (a'b')c' \text{ or } f' = a'(b'c')
@@ -380,7 +386,7 @@ Whereas computing from the right to left is first $jk1$ operations leaving a $j
* left to right is $njk + nk$ = $nk \cdot (j + 1)$.
* right to left is $jk + jn = j\cdot (k+n)$.
When $j=k$, say, we can compare and see the second is a factor less in terms of operations. This can be quite significant in higher dimensions, whereas the dimensions of calculus (where $n$ and $m$ are $3$ or less) it is not an issue.
When $j=k$, say, we can compare and see the second is a factor less in terms of operations. This can be quite significant in higher dimensions.
##### Example
@@ -405,7 +411,7 @@ Whereas the relationship is changed when the first matrix is skinny and the last
----
In calculus, we have $n$ and $m$ are $1$,$2$,or $3$. But that need not be the case, especially if differentiation is over a parameter space.
In calculus, we typically have $n$ and $m$ are $1$, $2$,or $3$. But that need not be the case, especially if differentiation is over a parameter space.
## Derivatives of matrix functions
@@ -546,7 +552,7 @@ Appropriate sizes for $A$, $B$, and $C$ are determined by the various products i
If $A$ is $m \times n$ and $B$ is $r \times s$, then since $BC$ is defined, $C$ has $s$ rows, and since $CA^T$ is defined, $C$ must have $n$ columns, as $A^T$ is $n \times m$, so $C$ must be $s\times n$. Checking this is correct on the other side, $A \times B$ would be size $mr \times ns$ and $\vec{C}$ would be size $sn$, so that product works, size wise.
The referred to notes have an explanation for this formula, but we confirm with an example with $m=n=2$ and $r=s=3$:
The referred to notes have an explanation for this formula, but we only confirm it with an example using $m=n=2$ and $r=s=3$:
```{julia}
@syms A[1:2, 1:2]::real B[1:3, 1:3]::real C[1:3, 1:2]::real

5
quarto/integral_vector_calculus/Project.toml
View File

@@ -4,9 +4,14 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
ImplicitIntegration = "bc256489-3a69-4a66-afc4-127cc87e6182"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
PlotlyKaleido = "f2990250-8cf9-495f-b13a-cce12b45703c"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
QuizQuestions = "612c44de-1021-4a21-84fb-7261cf5eb2d4"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
TextWrap = "b718987f-49a8-5099-9789-dcd902bef87d"

3
quarto/integrals/make_pdf.jl
View File

@@ -15,7 +15,8 @@ files = (
"center_of_mass",
"volumes_slice",
"arc_length",
"surface_area",
"surface_area",
"orthogonal_polynomials",
"twelve-qs",
)

724
quarto/integrals/orthogonal_polynomials.qmd Normal file
View File

@@ -0,0 +1,724 @@
# Orthogonal polynomials
{{< include ../_common_code.qmd >}}
This section uses these add-on packages:
```{julia}
using SymPy
using QuadGK
using Roots
using ForwardDiff: derivative
```
This section takes a detour to give some background on why the underlying method of `quadgk` is more efficient than those of Riemann sums. Orthogonal polynomials play a key role. There are many families of such polynomials. We highlight two.
## Inner product
Define an operation between two integrable, real-valued functions $f(x)$ and $g(x)$ by:
$$
\langle f, g \rangle = \int_{-1}^1 f(x)g(x) dx
$$
The properties of the integral mean this operation satisfies these three main properties:
* symmetry: $\langle f, g \rangle = \langle g,f \rangle$
* positive definiteness: $\langle f, f \rangle > 0$ *unless* $f(x)=0$.
* linearity: if $a$ and $b$ are scalars, then $\langle af + bg, h \rangle = a\langle f, h \rangle + b \langle g, h \rangle$.
The set of integrable functions forms a *vector space*, which simply means two such functions can be added to yield another integrable function and an integrable function times a scalar is still an integrable function. Many different collections of objects form a vector space. In particular, other sets of functions form a vector space, for example the collection of polynomials of degree $n$ or less or just the set of all polynomials.
For a vector space, an operation like the above satisfying these three properties is called an *inner product*; the combination of an inner product and a vector space is called an *inner product space*. In the following, we assume $f$ and $g$ are from a vector space with a real-valued inner product.
Inner products introduce a sense of size through a *norm*:
$\lVert f \rVert = \sqrt{\langle f, f\rangle }$.
Norms satisfy two main properties:
* scalar: $\lVert af \rVert = |a|\lVert f\rVert$
* triangle inequality: $\lVert f + g \rVert \leq \lvert f \rVert + \lVert g \rVert$
Two elements of an inner product space, $f$ and $g$, are *orthogonal* if $\langle f, g \rangle = 0$. This is a generalization of perpendicular. The Pythagorean theorem for orthogonal elements holds: $\lVert f\rVert^2 + \lVert g\rVert^2 = \lVert f+g\rVert^2$.
As we assume a real-valued inner product, the angle between two elements can be defined by:
$$
\angle(f,g) = \cos^{-1}\left(\frac{\langle f, g\rangle}{\lVert f \rVert \lVert g \rVert}\right).
$$
This says, the angle between two orthogonal elements is $90$ degrees (in some orientation)
The Cauchy-Schwarz inequality, $|\langle f, g \rangle| \leq \lVert f \rVert \lVert g\rVert$, for an inner product space, ensures the argument to $\cos^{-1}$ is between $-1$ and $1$.
These properties generalize two-dimensional vectors, with components $\langle x, y\rangle$. Recall, these can be visualized by placing a tail at the origin and a tip at the point $(x,y)$. Such vectors can be added by placing the tail of one at the tip of the other and using the vector from the other tail to the other tip.
With this, we have a vector anchored at the origin can be viewed as a line segment with slope $y/x$ (rise over run). A perpendicular line segment would have slope $-x/y$ (the negative reciprocal) which would be associated with the vector $\langle y, -x \rangle$. The dot product is just the sum of the multiplied components, or for these two vectors $x\cdot y + y\cdot (-x)$, which is $0$, as the line segments are perpendicular (orthogonal).
Consider now two vectors, say $f$, $g$. We can make a new vector that is orthogonal to $f$ by combining $g$ with a piece of $f$. But what piece?
Consider this
$$
\begin{align*}
\langle f, g - \frac{\langle f,g\rangle}{\langle f, f\rangle} f \rangle
&= \langle f, g \rangle - \langle f, \frac{\langle f,g\rangle}{\langle f, f\rangle} f \rangle \\
&= \langle f, g \rangle - \frac{\langle f,g\rangle}{\langle f, f\rangle}\langle f,f \rangle \\
&= \langle f, g \rangle - \langle f, g \rangle = 0
\end{align*}
$$
Define
$$
proj_f(g) = \frac{\langle f,g\rangle }{\langle f, f\rangle} f,
$$
then we have $u_1 = f$ and $u_2 = g-proj_f(g)$, $u_1$ and $u_2$ are orthogonal.
A similar calculation shows if $h$ is added to the set of elements, then
$u_3 = h - proj_{u_1}(h) - proj_{u_2}(h)$ will be orthogonal to $u_1$ and $u_2$. etc.
This process, called the [Gram-Schmidt](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process) process, can turn any set of vectors into a set of orthogonal vectors, assuming they are all non zero and no non-trivial linear combination makes them zero.
## Legendre
Consider now polynomials of degree $n$ or less with the normalization that $p(1) = 1$. We begin with two such polynomials: $u_0(x) = 1$ and $u_1(x) = x$.
These are orthogonal with respect to $\int_{-1}^1 f(x) g(x) dx$, as
$$
\int_{-1}^1 u_0(x) u_1(x) dx =
\int_{-1}^1 1 \cdot x dx =
x^2 \mid_{-1}^1 = 1^2 - (-1)^2 = 0.
$$
Now consider a quadratic polynomial, $u_2(x) = ax^2 + bx + c$, we want a polynomial which is orthogonal to $u_0$ and $u_1$ with the extra condition that $u_2(1) = c =1$ (or $c=1$.). We can do this using Gram-Schmidt as above, or as here through a system of two equations:
```{julia}
@syms a b c d x
u0 = 1
u1 = x
u2 = a*x^2 + b*x + c
eqs = (integrate(u0 * u2, (x, -1, 1)) ~ 0,
integrate(u1 * u2, (x, -1, 1)) ~ 0)
sols = solve(eqs, (a, b, c)) # b => 0, a => -3c
u2 = u2(sols...)
u2 = simplify(u2 / u2(x=>1)) # make u2(1) = 1 and fix c
```
The quadratic polynomial has $3$ unknowns and the orthgonality conditions give two equations. Solving these equations leaves one unknown (`c`). But the normalization condition (that $u_i(1) = 1$) allows `c` to be simplified out.
We can do this again with $u_3$:
```{julia}
u3 = a*x^3 + b*x^2 + c*x + d
eqs = (integrate(u0 * u3, (x, -1, 1)) ~ 0,
integrate(u1 * u3, (x, -1, 1)) ~ 0,
integrate(u2 * u3, (x, -1, 1)) ~ 0)
sols = solve(eqs, (a, b, c, d)) # a => -5c/3, b=>0, d=>0
u3 = u3(sols...)
u3 = simplify(u3/u3(x=>1)) # make u3(1) = 1
```
In theory, this can be continued up until any $n$. The resulting
polynomials are called the
[Legendre](https://en.wikipedia.org/wiki/Legendre_polynomials)
polynomials.
Rather than continue this, we develop easier means to generate these polynomials.
## General weight function
Let $w(x)$ be some non-negative function and consider the new inner product between two polynomials:
$$
\langle p, q\rangle = \int_I p(x) q(x) w(x) dx
$$
where $I$ is an interval and $w(x)$ is called a weight function. In the above discussion $I=[-1,1]$ and $w(x) = 1$.
Suppose we have *orthogonal* polynomials $p_i(x)$, $i=0,1, \dots, n$, where $p_i$ is a polynomial of degree $i$ ($p_i(x) = k_i x^i + \cdots$, where $k_i \neq 0$), and
$$
\langle p_m, p_n \rangle =
\int_I p_m(x) p_n(x) w(x) dx =
\begin{cases}
0 & m \neq n\\
h_m & m = n
\end{cases}
$$
Unique elements can be defined by specifying some additional property. For Legendre, it was $p_n(1)=1$, for other orthogonal families this may be specified by having leading coefficient of $1$ (monic), or a norm of $1$ (orthonormal), etc.
The above is the *absolutely continuous* case, generalizations of the integral allow this to be more general.
Orthogonality can be extended: If $q(x)$ is any polynomial of degree $m < n$, then
$\langle q, p_n \rangle = \int_I q(x) p_n(x) w(x) dx = 0$. (See the questions for more detail.)
Some names used for the characterizing constants are:
* $p_n(x) = k_n x^n + \cdots$ ($k_n$ is the leading term)
* $h_n = \langle p_n, p_n\rangle$
### Three-term reccurence
Orthogonal polynomials, as defined above through a weight function, satisfy a *three-term recurrence*:
$$
p_{n+1}(x) = (A_n x + B_n) p_n(x) - C_n p_{n-1}(x),
$$
where $n \geq 0$ and $p_{n-1}(x) = 0$.
(With this and knowledge of $A_n$, $B_n$, and $C_n$, the polynomials can be recursively generated from just specifying a value for the constant $p_0(x)$.
First, $p_{n+1}$ has leading term $k_{n+1}x^{n+1}$. Looking on the right hand side for the coefficient of $x^{n+1}$ we find $A_n k_n$, so $A_n = k_{n+1}/k_n$.
Next, we look at $u(x) = p_{n+1}(x) - A_n x p_n(x)$, a polynomial of degree $n$ or less.
As this has degree $n$ or less, it can be expressed in terms of $p_0, p_1, \dots, p_n$. Write it as $u(x) = \sum_{j=0}^n d_j p_j(x)$. Now, take any $m < n-1$ and consider $p_m$. We consider the inner product of $u$ and $p_m$ two ways:
$$
\begin{align*}
\int_I p_m(x) u(x) w(x) dx &=
\int_I p_m(x) \sum_{j=0}^n p_j(x) w(x) dx \\
&= \int_I p_m(x) \left(p_m(x) + \textcolor{red}{\sum_{j=0, j\neq m}^{n} p_j(x)}\right) w(x) dx \\
&= \int_I p_m(x) p_m(x) w(x) dx = h_m
\end{align*}
$$
As well
$$
\begin{align*}
\int_I p_m(x) u(x) w(x) dx
&= \int_I p_m(x) (p_{n+1}(x) - A_n x p_n(x)) w(x) dx \\
&= \int_I p_m(x) \textcolor{red}{p_{n+1}(x)} w(x) dx - \int_I p_m(x) A_n x p_n(x) w(x) dx\\
&= 0 - A_n \int_I (\textcolor{red}{x p_m(x)}) p_n(x) w(x) dx\\
&= 0
\end{align*}
$$
The last integral being $0$ as $xp_m(x)$ has degree $n-1$ or less and hence is orthogonal to $p_n$.
That is $p_{n+1} - A_n x p_n(x) = d_n p_n(x) + d_{n-1} p_{n-1}(x)$. Setting $B_n=d_n$ and $C_{n-1} = -d_{n-1}$ shows the three-term recurrence applies.
#### Example: Legendre polynomials
With this notation, the Legendre polynomials have:
$$
\begin{align*}
w(x) &= 1\\
I &= [-1,1]\\
A_n &= \frac{2n+1}{n+1}\\
B_n &= 0\\
C_n & = \frac{n}{n+1}\\
k_{n+1} &= \frac{2n+1}{n+1}k_n - \frac{n}{n-1}k_{n-1}, k_1=k_0=1\\
h_n &= \frac{2}{2n+1}
\end{align*}
$$
#### Favard theorem
In an efficient review of the subject, [Koornwinder](https://arxiv.org/pdf/1303.2825) states conditions on the recurrence that ensure that if a $n$-th degree polynomials $p_n$ satisfy a three-term recurrence, then there is an associated weight function (suitably generalized). The conditions use this form of a three-term recurrence:
$$
\begin{align*}
xp_n(x) &= a_n p_{n+1}(x) + b_n p_n(x) + c_n p_{n-1}(x),\quad (n > 0)\\
xp_0(x) &= a_0 p_1(x) + b_0 p_0(x)
\end{align*}
$$
where the constants are real and $a_n c_{n+1} > 0$. These force $a_n = k_n/k_{n+1}$ and $c_n/h_{n+1} = a_n/h_n$
#### Clenshaw algorithm
When introducing polynomials, the synthetic division algorithm was given to compute $p(x) / (x-r)$. This same algorithm also computed $p(r)$ efficiently and is called Horner's method. The `evalpoly` method in `Julia`'s base implements this.
For a set of polynomials $p_0(x), p_1(x), \dots, p_n(x)$ satisfying a three-term recurrence $p_{n+1}(x) = (A_n x + B_n) p_n(x) - C_n p_{n-1}(x)$, the Clenshaw algorithm gives an efficient means to compute an expression of a linear combination of the polynomials, $q(x) = a_0 p_0(x) + a_1 p_1(x) + \cdots + a_n p_n(x)$.
The [method](https://en.wikipedia.org/wiki/Clenshaw_algorithm) uses a reverse recurrence formula starting with
$$
b_{n+1}(x) = b_{n+2}(x) = 0
$$
and then computing for $k = n, n-1, \dots, 1$
$$
b_k(x) = a_k + (A_k x + B_k) b_{k+1}(x) - C_k b_{k+2}(x)
$$
Finally finishing by computing $a_0 p_0(x) + b_1 p_1(x) - C(1) p_0(x) b_2$.
For example, with the Legendre polynomials, we have
```{julia}
A(n) = (2n+1)//(n+1)
B(n) = 0
C(n) = n // (n+1)
```
Say we want to compute $a_0 u_0(x) + a_1 u_1(x) + a_2 u_2(x) + a_3 u_3(x) + a_4 u_4(x)$. The necessary inputs are the coefficients, the value of $x$, and polynomials $p_0$ and $p_1$.
```{julia}
function clenshaw(x, as, p0, p1)
n = length(as) - 1
bn1, bn2 = 0, 0
a(k) = as[k + 1] # offset
for k in n:-1:1
bn1, bn2 = a(k) + (A(k) * x + B(k)) * bn1 - C(k+1) * bn2, bn1
end
b1, b2 = bn1, bn2
p0(x) * a(0) + p1(x) * b1 - C(1) * p0(x) * b2
end
```
This function can be purposed to generate additional Legendre polynomials. For example, to compute $u_4$ we pass in a symbolic value for $x$ and mask out all by $a_4$ in our coefficients:
```{julia}
p₀(x) = 1
p₁(x) = x # Legendre
@syms x
clenshaw(x, (0,0,0,0,1), p₀, p₁) |> expand |> simplify
```
:::{.callout-note}
### Differential equations approach
A different description of families of orthogonal polynomials is that they satisfy a differential equation of the type
$$
\sigma(x) y''(x) + \tau(x) y'(x) + \lambda_n y(x) = 0,
$$
where $\sigma(x) = ax^2 + bx + c$, $\tau(x) = dx + e$, and $\lambda_n = -(a\cdot n(n-1) + dn)$.
With this parameterization, values for $A_n$, $B_n$, and $C_n$ can be given in terms of the leading coefficient, $k_n$ (cf. [Koepf and Schmersau](https://arxiv.org/pdf/math/9612224)):
$$
\begin{align*}
A_n &= \frac{k_{n+1}}{k_n}\\
B_n &= \frac{k_{n+1}}{k_n} \frac{2bn(a(n-1)+d) + e(d-2a)}{(2a(n-1) + d)(2an+d)}\\
C_n &= \frac{k_{n+1}}{k_{n-1}}
\frac{n(a(n-1) + d)(a(n-2)+d)(n(an+d))(4ac-b^2)+ae^2+cd^2-bde}{
(a(n-1)+d)(a(2n-1)+d)(a(2n-3)+d)(2a(n-1)+d)^2}
\end{align*}
$$
There are other relations between derivatives and the orthogonal polynomials. For example, another three-term recurrence is:
$$
\sigma(x) p_n'(x) = (\alpha_n x + \beta_n)p_n(x) + \gamma_n p_{n-1}(x)
$$
The same reference has formulas for $\alpha$, $\beta$, and $\gamma$ in terms of $a,b,c,d$, and $e$ along with many others.
:::
## Chebyshev
The Chebyshev polynomials (of the first kind) satisfy the three-term recurrence
$$
T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)
$$
with $T_0(x)= 1$ and $T_1(x)=x$.
These polynomials have domain $(-1,1)$ and weight function $(1-x^2)^{-1/2}$.
(The Chebyshev polynomials of the second kind satisfy the same three-term recurrence but have $U_0(x)=1$ and $U_1(x)=2x$.)
These polynomials are related to trigonometry through
$$
T_n(\cos(\theta)) = \cos(n\theta)
$$
This characterization makes it easy to find the zeros of the
polynomial $T_n$, as they happen when $\cos(n\theta)$ is $0$, or when
$n\theta = \pi/2 + k\pi$ for $k$ in $0$ to $n-1$. Solving for $\theta$
and taking the cosine, we get the zeros of the $n$th degree polynomial
$T_n$ are $\cos(\pi(k + 1/2)/n)$ for $k$ in $0$ to $n-1$.
These evenly spaced angles lead to roots more concentrated at the edges of the interval $(-1,1)$.
##### Example
Chebyshev polynomials have a minimal property that makes them fundamental for use with interpolation.
Define the *infinity* norm over $[-1,1]$ to be the maximum value of the absolute value of the function over these values.
Let $f(x) = 2^{-n+1}T_n(x)$ be a monic version of the Chebyshev polynomial.
> If $q(x)$ is any monic polynomial of degree $n$, then the infinity norm of $q(x)$ is greater than or equal to that of $f$.
Using the trigonometric representation of $T_n$, we have
* $f(x)$ has infinity norm of $2^{-n+1}$ and these maxima occur at $x=\cos((k\pi)/n)$, where $0 \leq k \leq n$. (There is a cosine curve with known peaks, oscillating between $-1$ and $1$.)
* $f(x) > 0$ at $x = \cos((2k\pi)/n)$ for $0 \leq 2k \leq n$
* $f(x) < 0$ at $x = \cos(((2k+1)\pi)/n)$ for $0 \leq 2k+1 \leq n$
Suppose $w(x)$ is a monic polynomial of degree $n$ and suppose it has smaller infinity norm. Consider $u(x) = f(x) - w(x)$. At the extreme points of $f(x)$, we must have $|f(x)| \geq |w(x)|$. But this means
* $u(x) > 0$ at $x = \cos((2k\pi)/n)$ for $0 \leq 2k \leq n$
* $u(x) < 0$ at $x = \cos(((2k+1)\pi)/n)$ for $0 \leq 2k+1 \leq n$
As $u$ is continuous, this means there are at least $n$ sign changes, hence $n$ or more zeros. But as both $f$ and $w$ are monic, $u$ is of degree $n-1$, at most. This is a contradiction unless $u(x)$ is the zero polynomial, which it can't be by assumption.
### Integration
Recall, a Riemann sum can be thought of in terms of weights, $w_i$ and nodes $x_i$ for which $\int_I f(x) dx \approx \sum_{i=0}^{n-1} w_i f(x_i)$.
For a right-Riemann sum with partition given by $a_0 < a_1 < \cdots < a_n$ the nodes are $x_i = a_i$ and the weights are $w_i = (a_i - a_{i-1})$ (or in the evenly spaced case, $w_i = (a_n - a_0)/n$.
More generally, this type of expression can represent integrals of the type $\int_I f(x) w(x) dx$, with $w(x)$ as in an inner product. Call such a sum a Gaussian quadrature.
We will see that the zeros of orthogonal polynomials have special properties as nodes.
> For orthogonal polynomials over the interval $I$ with weight function $w(x)$, each $p_n$ has $n$ distinct real zeros in $I$.
Suppose that $p_n$ had only $k<n$ sign changes at $x_1, x_2, \dots, x_k$. Then for some choice of $\delta$, $(-1)^\delta p(x) (x-x_1)(x-x_2)\cdots(x-x_k) \geq 0$. Since this is non zero, it must be that
$$
(-1)^\delta \int_I p(x) \left( (x-x_1)(x-x_2)\cdots(x-x_k)\right) w(x) dx > 0
$$
But, the product is of degree $k < n$, so by orthogonality must be $0$. Hence, it can't be that $k < n$, so there must be $n$ sign changes in $I$ by $p_n$. Each corresponds to a zero, as $p_n$ is continuous.
This next statement says that using the zeros of $p_n$ for the nodes of Gaussian quadrature and appropriate weights that the quadrature is exact for higher degree polynomials.
> For a fixed $n$, suppose $p_0, p_1, \dots, p_n$ are orthogonal polynomials over $I$ with weight function $w(x)$. If the zeros of $p_n$ are the nodes $x_i$, then there exists $n$ weights so that the any polynomial of degree $2n-1$ or less, the Gaussian quadrature is exact.
That is if $q(x)$ is a polynomial with degree $2n-1$ or less, we have for some choice of $w_i$:
$$
\int_I q(x) w(x) dx = \sum_{i=1}^n w_i q(x_i)
$$
To compare, recall, Riemann sums ($1$-node) were exact for constant functions (degree $0$), the trapezoid rule ($2$-nodes) is exact for linear polynomials (degree $1$), and Simpson's rule ($3$ nodes) are exact for cubic polynomials (degree $3$).
We follow [Wikipedia](https://en.wikipedia.org/wiki/Gaussian_quadrature#Fundamental_theorem) to see this key fact.
Take $h(x)$ of degree $2n-1$ or less. Then by polynomial long division, there are polynomials $q(x)$ and $r(x)$ where
$$
h(x) = q(x) p_n(x) + r(x)
$$
and the degree of $r(x)$ is less than $n-1$, the degree of $p_n(x)$. Further, the degree of $q(x)$ is also less than $n-1$, as were it more, then the degree of $q(x)p_n(x)$ would be more than $n-1+n$ or $2n-1$. Let's note that if $x_i$ is a zero of $p_n(x)$ that $h(x_i)= r(x_i)$.
So
$$
\begin{align*}
\int_I h(x) w(x) dx &= \int_I \textcolor{red}{q(x)} p_n(x) w(x) dx + \int_I r(x) w(x)dx\\
&= 0 + \int r(x) w(x) dx.
\end{align*}
$$
Now consider the polynomials made from the zeros of $p_n(x)$
$$
l_i(x) = \prod_{j \ne i} \frac{x - x_j}{x_i - x_j}
$$
These are called Lagrange interpolating polynomials and have the property that $l_i(x_i) = 1$ and $l_i(x_j) = 0$ if $i \neq j$.
These allow the expression of
$$
\begin{align*}
r(x) &= l_1(x)r(x_1) + l_2(x) r(x_2) + \cdots + l_n(x) r(x_n) \\
&= \sum_{i=1}^n l_i(x) r(x_i)
\end{align*}
$$
This isn't obviously true, but this expression agrees with an at-most degree $n-1$ polynomial ($r(x)$) at $n$ points hence it must be the same polynomial.)
With this representation, the integral becomes
$$
\begin{align*}
\int_I h(x) w(x) dx &= \int_I r(x) w(x)dx \\
&= \int_I \sum_{i=1}^n l_i(x) r(x_i) w(x) dx\\
&= \sum_{i=1}^n r(x_i) \int_I l_i(x) w(x) dx \\
&= \sum_{i=1}^n r(x_i) w_i\\
&= \sum_{i=1}^n w_i h(x_i)
\end{align*}
$$
That is there are weights, $w_i = \int_I l_i(x) w(x) dx$, for which the integration is exactly found by Gaussian quadrature using the roots of $p_n$ as the nodes.
The general formula for the weights can be written in terms of the polynomials $p_i = k_ix^i + \cdots$:
$$
\begin{align*}
w_i &= \int_I l_i(x) w(x) dx \\
&= \frac{k_n}{k_{n-1}}
\frac{\int_I p_{n-1}(x)^2 w(x) dx}{p'_n(x_i) p_{n-1}(x_i)}.
\end{align*}
$$
To see this, consider:
$$
\begin{align*}
\prod_{j \neq i} (x - x_j) &=
\frac{\prod_j (x-x_j)}{x-x_i} \\
&= \frac{1}{k_n}\frac{k_n \prod_j (x - x_j)}{x - x_i} \\
&= \frac{1}{k_n} \frac{p_n(x)}{x-x_i}\\
&= \frac{1}{k_n} \frac{p_n(x) - p_n(x_i)}{x-x_i}\\
&\rightarrow \frac{p'_n(x_i)}{k_n}, \text{ as } x \rightarrow x_i.
\end{align*}
$$
Thus
$$
\prod_{j \neq i} (x_i - x_j) = \frac{p'_n(x_i)}{k_n}.
$$
This gives
$$
\begin{align*}
w_i &= \int_i \frac{k_n \prod_j (x-x_j)}{p'_n(x_i)} w(x) dx\\
&= \frac{1}{p'_n(x_i)} \int_i \frac{p_n(x)}{x-x_i} w(x) dx
\end{align*}
$$
To work on the last term, a trick (see the questions for detail) can show that for any $k \leq n$ that
$$
\int_I \frac{x^k p_n(x)}{x - x_i} w(x) dx
= x_i^k \int_I \frac{p_n(x)}{x - x_i} w(x) dx
$$
Hence for any degree $n$ or less polynomial: we have
$$
q(x_i) \int_I \frac{p_n(x)}{x - x_i} w(x) dx =
\int_I \frac{q(x) p_n(x)}{x - x_i} w(x) dx
$$.
We will use this for $p_{n-1}$. First, as $x_i$ is a zero of $p_n(x)$ we have
$$
\frac{p_n(x)}{x-x_i} = k_n x^{n-1}+ r(x),
$$
where $r(x)$ has degree $n-2$ at most. This is due to $p_n$ being divided by a monic polynomial, hence leaving a degree $n-1$ polynomial with leading coefficient $k_n$.
But then
$$
\begin{align*}
w_i &= \frac{1}{p'_n(x_i)} \int_I \frac{p_n(x)}{x-x_i} w(x) dx \\
&= \frac{1}{p'_n(x_i)} \frac{1}{p_{n-1}(x_i)} \int_I \frac{p_{n-1}(x) p_n(x)}{x - x_i} w(x) dx\\
&= \frac{1}{p'_n(x_i)p_{n-1}(x_i)} \int_I p_{n-1}(x)
(k_n x^{n-1} + \textcolor{red}{r(x)}) w(x) dx\\
&= \frac{k_n}{p'_n(x_i)p_{n-1}(x_i)} \int_I p_{n-1}(x) x^{n-1} w(x) dx\\
&= \frac{k_n}{p'_n(x_i)p_{n-1}(x_i)} \int_I p_{n-1}(x)
\left(
\textcolor{red}{\left(x^{n-1} - \frac{p_{n-1}(x)}{k_{n-1}}\right) }
+ \frac{p_{n-1}(x)}{k_{n-1}}\right) w(x) dx\\
&= \frac{k_n}{p'_n(x_i)p_{n-1}(x_i)} \int_I p_{n-1}(x)\frac{p_{n-1}(x)}{k_{n-1}} w(x) dx\\
&= \frac{k_n}{k_{n-1}} \frac{1}{p'_n(x_i)p_{n-1}(x_i)} \int_I p_{n-1}(x)^2 w(x) dx.
\end{align*}
$$
### Examples of quadrature formula
The `QuadGK` package uses a modification to Gauss quadrature to estimate numeric integrals. Let's see how. Behind the scenes, `quadgk` calls `kronrod` to compute nodes and weights.
We have from earlier that
```{julia}
u₃(x) = x*(5x^2 - 3)/2
u₄(x) = 35x^4 / 8 - 15x^2 / 4 + 3/8
```
```{julia}
xs = find_zeros(u₄, -1, 1)
```
From this we can compute the weights from the derived general formula:
```{julia}
k₃, k₄ = 5/2, 35/8
w(x) = 1
I = first(quadgk(x -> u₃(x)^2 * w(x), -1, 1))
ws = [k₄/k₃ * 1/(derivative(u₄,xᵢ) * u₃(xᵢ)) * I for xᵢ ∈ xs]
(xs, ws)
```
We compare now to the values returned by `kronrod` in `QuadGK`
```{julia}
kxs, kwts, wts = kronrod(4, -1, 1)
[ws wts xs kxs[2:2:end]]
```
(The `kronrod` function computes $2n-1$ nodes and weights. The Gauss-Legendre nodes are $n$ of those, and extracted by taking the 2nd, 4th, etc.)
To compare integrations of some smooth function we have
```{julia}
u(x) = exp(x)
GL = sum(wᵢ * u(xᵢ) for (xᵢ, wᵢ) ∈ zip(xs, ws))
KL = sum(wᵢ * u(xᵢ) for (xᵢ, wᵢ) ∈ zip(kxs, kwts))
QL, esterror = quadgk(u, -1, 1)
(; GL, KL, QL, esterror)
```
The first two are expected to not be as accurate, as they utilize a fixed number of nodes.
## Questions
###### Question
Let $p_i$ for $i$ in $0$ to $n$ be polynomials of degree $i$. It is true that for any polynomial $q(x)$ of degree $k \leq n$ that there is a linear combination such that $q(x) = a_0 p_0(x) + \cdots + a_k p_k(x)$.
First it is enough to do this for a monic polynomial $x^k$, why?
```{julia}
#| echo: false
choices = [raw"If you can do it for each $x^i$ then if $q(x) = b_0 + b_1x + b_2x^2 + \cdots + b_k x^k$ we just multiply the coefficients for each $x^i$ by $b_i$ and add.",
raw"It isn't true"]
radioq(choices, 1)
```
Suppose $p_0 = k_0$ and $p_1 = k_1x + a$. How would you make $x=x^1$?
```{julia}
#| echo: false
choices = [raw"$(p1 - (a/k_0) p_0)/k_1$",
raw"$p1 - p0$"]
radioq(choices, 1)
```
Let $p_i = k_i x^i + \cdots$ ($k_i$ is the leading term)
To reduce $p_3 = k_3x^3 + a_2x^2 + a_1x^1 + a_0$ to $k_3x^3$ we could try:
* form $q_3 = p_3 - p_2 (a_2/k_2)$. As $p_2$ is degree $2$, this leaves $k_3x^3$ alone, but it
```{julia}
#| echo: false
choices = [raw"It leaves $0$ as the coefficient of $x^2$",
raw"It leaves all the other terms as $0$"]
radioq(choices, 1)
```
* We then use $p_1$ times some multiple $a/k_1$ to remove the $x$ term
* we then use $p_0$ times some multiple $a/k_0$ to remove the constant term
Would this strategy work to reduce $p_n$ to $k_n x^n$?
```{julia}
#| echo: false
radioq(["Yes", "No"], 1)
```
###### Question
Suppose $p(x)$ and $q(x)$ are polynomials of degree $n$ and there are $n+1$ points for which $p(x_i) = q(x_i)$.
First, is it true or false that a polynomial of degree $n$ has *at most* n zeros?
```{julia}
#| echo: false
radioq(["true, unless it is the zero polynomial", "false"], 1)
```
What is the degree of $h(x) = p(x) - q(x)$?
```{julia}
#| echo: false
radioq([raw"At least $n+1$", raw"At most $n$"], 2)
```
At least how many zeros does the polynomial $h(x)$ have?
```{julia}
#| echo: false
radioq([raw"At least $n+1$", raw"At most $n$"], 1)
```
Is $p(x) = q(x)$ with these assumptions?
```{julia}
#| echo: false
radioq(["yes", "no"], 1)
```
###### Question
We wish to show that for any $k \leq n$ that
$$
\int_I \frac{x^k p_n(x)}{x - x_i} w(x) dx
= x_i^k \int_I \frac{p_n(x)}{x - x_i} w(x) dx
$$
We have for $u=x/x_i$ that
$$
\frac{1}{x - x_i} = \frac{1 - u^k}{x - x_i} + \frac{u^k}{x - x_i}
$$
But the first term, $(1-u^k)/(x-x_i)$ is a polynomial of degree $k-1$. Why?
```{julia}
#| echo: false
choices = [raw"""
Because we can express this as $x_i^k - x^k$ which factors as $(x_i - x) \cdot u(x)$ where $u(x)$ has degree $k-1$, at most.
""",
raw"""
It isn't true, it clearly has degree $k$
"""]
radioq(choices, 1)
```
This gives if $k \leq n$ and with $u=x/x_i$:
$$
\begin{align*}
\int_I \frac{p_n(x)}{x - x_i} w(x) dx
&= \int_I p_n(x) \left( \textcolor{red}{\frac{1 - u^k}{x - x_i}} + \frac{u^k}{x - x_i} \right) w(x) dx\\
&= \int_I p_n(x) \frac{\frac{x^k}{x_i^k}}{x - x_i} w(x) dx\\
&= \frac{1}{x_i^k} \int_I \frac{x^k p_n(x)}{x - x_i} w(x) dx
\end{align*}
$$

12
quarto/misc/Project.toml Normal file
View File

@@ -0,0 +1,12 @@
[deps]
CalculusWithJulia = "a2e0e22d-7d4c-5312-9169-8b992201a882"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
QuizQuestions = "612c44de-1021-4a21-84fb-7261cf5eb2d4"
SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
TextWrap = "b718987f-49a8-5099-9789-dcd902bef87d"