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

WIP; better figures

Browse Source
This commit is contained in:
jverzani
2025-07-02 14:05:09 -04:00
parent 5013211954
commit c38a7c9f1d
5 changed files with 380 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
quarto/adjust_plotly.jl
View File

@@ -9,7 +9,6 @@
function _add_plotly(f)
#lineno = 117
r = readlines(f)
inserted = false
open(f, "w") do io

31
quarto/integrals/arc_length.qmd
View File

@@ -76,21 +76,34 @@ To see why, any partition of the interval $[a,b]$ by $a = t_0 < t_1 < \cdots < t
## {{{arclength_graph}}}
gr()
function make_arclength_graph(n)
x(t) = cos(t)/t
y(t) = sin(t)/t
a, b = 1, 4pi
ns = [10,15,20, 30, 50]
empty_style = (xaxis=([], false),
yaxis=([], false),
framestyle=:origin,
legend=false)
ns = [10,15,20, 30, 50]
plot(; empty_style..., aspect_ratio=:equal, size=fig_size)
title!("Approximate arc length with $(ns[n]) points")
g(t) = cos(t)/t
f(t) = sin(t)/t
ts = range(a, b, 250)
plot!(x.(ts), y.(ts); line=(:black,2))
pttn = range(a, b, ns[n])
plot!(x.(pttn), y.(pttn); line=(:red, 2))
ts = range(1, stop=4pi, length=200)
tis = range(1, stop=4pi, length=ns[n])
ts = range(0, 2pi, 100)
p = plot(g, f, 1, 4pi, legend=false, size=fig_size,
title="Approximate arc length with $(ns[n]) points")
plot!(p, map(g, tis), map(f, tis), color=:orange)
λ = 0.005
cs = [λ .* xys for xys ∈ sincos.(ts)]
p
for v ∈ zip(x.(pttn), y.(pttn))
S = Shape([v .+ xy for xy in cs])
plot!(S; fill=(:white,), line=(:black,2))
end
current()
end
n = 5

69
quarto/integrals/area.qmd
View File

@@ -190,6 +190,40 @@ Clearly for a given partition and choice of $c_i$, the above can be computed. Ea
#| hold: true
#| echo: false
gr()
function left_riemann(n)
empty_style = (xaxis=([], false),
yaxis=([], false),
framestyle=:origin,
legend=false)
axis_style = (arrow=true, side=:head, line=(:gray, 1))
rectangle = (x, y, w, h) -> Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
f = x -> -(x+1/2)*(x-1)*(x-3) + 1
a, b= 1, 3
plot(; empty_style...)
plot!(f, a, b; line=(:black, 3))
plot!([a-.25, b+.25], [0,0]; axis_style...)
plot!([a-.1, a-.1], [-.25, .5 + f(a/2 +b/2)]; axis_style...)
Δ = (b-a)/n
for i ∈ 0:n-1
xᵢ = a + i*Δ
plot!(rectangle(xᵢ, 0, Δ, f(xᵢ)), opacity=0.5, color=:red)
end
area = round(sum(f(a + i*Δ)*Δ for i ∈ 0:n-1), digits=3)
annotate!([
(a, 0, text(L"a", :top)),
(b, 0, text(L"b", :top)),
(a, f(a/2+b/2), text("\$L_{$n} = $area\$", :left))
])
current()
end
#=
rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
function ₙ(j)
a = ("₋","","","₀","₁","₂","₃","₄","₅","₆","₇","₈","₉")
@@ -210,6 +244,7 @@ function left_riemann(n)
title!("L$(ₙ(n)) = $a")
p
end
=#
anim = @animate for i ∈ (2,4,8,16,32,64)
left_riemann(i)
@@ -631,19 +666,19 @@ Before continuing, we define a function to compute Riemann sums for us with an
```{julia}
#| eval: false
riemann(f, a, b, n; method="right") = riemann(f, range(a,b,n+1); method=method)
function riemann(f, xs; method="right")
Ms = (left = (f,a,b) -> f(a),
right= (f,a,b) -> f(b),
Ms = (left = (f,a,b) -> f(a),
right = (f,a,b) -> f(b),
trapezoid = (f,a,b) -> (f(a) + f(b))/2,
simpsons = (f,a,b) -> (c = a/2 + b/2; (1/6) * (f(a) + 4*f(c) + f(b))),
simpsons = (f,a,b) -> (c = a/2 + b/2; (1/6) * (f(a) + 4*f(c) + f(b)))
)
_riemann(Ms[Symbol(method)], f, xs)
end
function _riemann(M, f, xs)
M = Ms[Symbol(method)}
xs′ = zip(xs[1:end-1], xs[2:end])
sum(M(f, a, b) * (b-a) for (a,b) ∈ xs′)
end
riemann(f, a, b, n; method="right") =
riemann(f, range(a,b,n+1); method)
```
(This function is defined in `CalculusWithJulia` and need not be copied over if that package is loaded.)
@@ -699,7 +734,7 @@ Consider a function $g(x)$ defined through its piecewise linear graph:
```{julia}
#| echo: false
g(x) = abs(x) > 2 ? 1.0 : abs(x) - 1.0
plot(g, -3,3)
plot(g, -3,3; legend=false)
plot!(zero)
```
@@ -710,6 +745,24 @@ plot!(zero)
We could add the signed area over $[0,1]$ to the above, but instead see a square of area $1$, a triangle with area $1/2$ and a triangle with signed area $-1$. The total is then $1/2$.
This figure may make the above decomposition more clear:
```{julia}
#| echo: false
let
g(x) = abs(x) > 2 ? 1.0 : abs(x) - 1.0
xs = [ -3, -2, -1, 1, 2, 3]
plot(; legend=false, aspect_ratio=:equal)
plot!(Shape([-3,-2,-2,-3], [0,0,1,1]); fill=(:gray,))
plot!(Shape([-2,-1,-1,-2], [0,0,0,1]); fill=(:gray10,))
plot!(Shape([-1,0,1], [0,-1,0]); fill=(:gray90,))
plot!(Shape([1,2,2,1], [0,1,0,0]); fill=(:gray10,))
plot!(Shape([2,3,3,2], [0,0,1,1]); fill=(:gray,))
end
```
* Compute $\int_{-3}^{3} g(x) dx$:

40
quarto/integrals/twelve-qs.qmd
View File

@@ -6,7 +6,12 @@ This section uses these packages:
using SymPy
using Plots
using Roots
plotly()
```
```{julia}
#| echo: false
using LaTeXStrings
gr();
```
----
@@ -15,21 +20,31 @@ In the March 2003 issue of the College Mathematics Journal, Leon M Hall posed 12
```{julia}
#| echo: false
f(x) = x^2
fp(x) = 2x
a₀ = 7/8
q₀ = -a₀ - 1/(2a₀)
f(x) = x^2
fp(x) = 2x
tangent(x) = f(a₀) + fp(a₀) * (x - a₀)
normal(x) = f(a₀) - (1 / fp(a₀)) * (x - a₀)
function make_plot(a₀=7/8, q₀=-a₀ - 1/2a₀)
plt = plot(; xlim=(-2,2), ylim=(-1, (1.5)^2),
xticks=nothing, yticks=nothing,
aspect_ratio=:equal, border=:none, legend=false)
function make_plot()
empty_style = (xaxis=([], false),
yaxis=([], false),
framestyle=:origin,
legend=false)
axis_style = (arrow=true, side=:head, line=(:gray, 1))
tangent(x) = f(a₀) + fp(a₀) * (x - a₀)
normal(x) = f(a₀) - (1 / fp(a₀)) * (x - a₀)
plt = plot(; empty_style...,
xlims=(-2,2), ylims=(-1, (1.5)^2))
f(x) = x^2
fp(x) = 2x
plot!(f, -1.5, 1.5, line=(1, :royalblue))
plot!(zero, line=(1, :black))
plot!(f, -1.5, 1.5, line=(2, :black))
plot!([-1.6, 1.6], [0,0]; axis_style...)
tl = x -> f(a₀) + fp(a₀) * (x-a₀)
nl = x -> f(a₀) - 1/(fp(a₀)) * (x-a₀)
@@ -40,9 +55,10 @@ function make_plot(a₀=7/8, q₀=-a₀ - 1/2a₀)
# add in right triangle
scatter!([a₀, q₀], f.([a₀, q₀]), markersize=5)
Δ = 0.01
annotate!([(a₀ + Δ, nl(a₀+Δ), "P", :bottom),
(q₀ - Δ, nl(q₀-Δ), "Q", :top)])
plt
annotate!([(a₀ + Δ, nl(a₀+Δ), text(L"P", :top)),
(q₀ - Δ, nl(q₀-Δ), text(L"Q", :bottom, :left))
])
current()
end
make_plot()
```

290
quarto/integrals/volumes_slice.qmd
View File

@@ -680,9 +680,93 @@ For a general cone, we use this [definition](http://en.wikipedia.org/wiki/Cone):
Let $h$ be the distance from the apex to the base. Consider cones with the property that all planes parallel to the base intersect the cone with the same shape, though perhaps a different scale. This figure shows an example, with the rays coming from the apex defining the volume.
::: {#fig-generic-cone}
```{julia}
#| echo: false
plt = let
gr()
rad(t) = 3/2 - t
viewp = [2,-1/1.5,1/2+.2]
empty_style = (xaxis=([], false),
yaxis=([], false),
framestyle=:origin,
legend=false)
axis_style = (arrow=true, side=:head, line=(:gray, 1))
unitize(x) = x / norm(x)
"""Orthogonal projection along the vector viewp"""
function make_Pmat(viewp)
a = unitize( [-viewp[2], viewp[1], 0] )
b = unitize( [-viewp[3]*viewp[1],
-viewp[3]*viewp[2],
viewp[1]^2 + viewp[2]^2] )
collect(zip(a,b))
end
#linear projection of R^3 onto R^2
function proj(X, viewp)
Pmat = make_Pmat(viewp)
x=sum([Pmat[i][1]*X[i] for i in 1:3])
y=sum([Pmat[i][2]*X[i] for i in 1:3])
(x, y) # a point
end
proj(X) = proj(X, viewp)
drawdiscF(t) = Shape([psurf(t, 2*i*pi/100) for i in 1:101])
# discrete determinant of Jacobian; area multiplier?
function jac(X, u, v)
ϵ = 0.000001
A = map((p,q) -> (p-q)/ϵ, X(u+ϵ/2, v), X(u-ϵ/2, v))
B = map((p,q) -> (p-q)/ϵ, X(u, v+ϵ/2), X(u, v-ϵ/2))
return A[1]*B[2]-A[2]*B[1]
end
# our surface
R, r, rho = 1, 1/4, 1/4
f(t) = (R-r) * cos(t) + rho * cos((R-r)/r * t)
g(t) = (R-r) * sin(t) - rho * sin((R-r)/r * t)
surf(t, θ) = (rad(t)*f(θ), rad(t)*g(θ), t)
psurf(t,θ) = proj(surf(t,θ))
plot(; empty_style..., aspect_ratio=:equal)
for (i,t) in enumerate(range(0, 3/2, 30))
plot!(drawdiscF(t); fill=(:gray,1), line=(:black,1))
end
θ = 0; plot!([psurf(0, θ), psurf(3/2, θ)]; line=(:black, 2))
θ = pi/2; plot!([psurf(0, θ), psurf(3/2, θ)]; line=(:black, 1))
θ = 3pi/2; plot!([psurf(0, θ), psurf(3/2, θ)]; line=(:black, 1))
current()
end
plt
```
```{julia}
#| echo: false
plotly()
nothing
```
A "cone" formed from the parameterized curve
$r(t) = \langle
(R-r) \cdot \cos(t) + \rho \cdot \cos((R-r)/r \cdot t),
(R-r) \cdot \sin(t) - \rho \cdot \sin((R-r)/r \cdot t)
\rangle$ with apex at the point $[0,0,3/2]$ and rays extending down through the origin following $3/2-z$.
:::
```{julia}
#| echo: false
#| eval: false
plt = let
h = 5
R, r, rho = 1, 1/4, 1/4
@@ -696,7 +780,7 @@ plt = let
end
current()
end
plot
plt
```
A right circular cone is one where this shape is a circle. This definition can be more general, as a square-based right pyramid is also such a cone. After possibly reorienting the cone in space so the base is at $u=0$ and the apex at $u=h$ the volume of the cone can be found from:
@@ -736,6 +820,11 @@ plt = let
r(h,z) = (x(h,z), y(h,z))
r1(h,z) = (2,0) .+ r(h,z)
empty_style = (xaxis=([], false),
yaxis=([], false),
framestyle=:origin,
legend=false)
Nh=30
heights = range(-1/2, 1/2, Nh)
h0=heights[Nh ÷ 2]
@@ -779,7 +868,8 @@ end
plt
```
This figure shows the volume of a figure being comprised of slices. A discrete approximation would be found by estimating the volume of each slice by the cross sectional area times a small $\Delta h$.
This figure shows the volume of a figure being comprised of slices. A discrete approximation would be found by estimating the volume of each slice by the cross sectional area times a small $\Delta h$. This leads to a formula
$V = \int_a^b A(h)dh$, where $A$ computes the cross sectional area.
(This figure was ported from @Angenent.)
:::
@@ -836,39 +926,197 @@ Illustration of Cavalieri's first principle. The discs from the left are moved a
:::
With the formula for the volume of solids based on cross sections, this is a trivial observation, as the functions giving the cross-sectional area are identical. Still, it can be surprising. Consider a sphere with an interior cylinder bored out of it. (The [Napkin](http://tinyurl.com/o237v83) ring problem.) The bore has height $h$ - for larger radius spheres this means very wide bores.
With the formula for the volume of solids based on cross sections, this is a trivial observation, as the functions giving the cross-sectional area are identical. Still, it can be surprising.
Consider a sphere with an interior cylinder bored out of it. (The [Napkin](http://tinyurl.com/o237v83) ring problem.) The bore has height $h$ - for larger radius spheres this means very wide bores.
::: {#fig-napkin-ring-1}
```{julia}
#| echo: false
plt = let
# Follow lead of # https://github.com/SigurdAngenent/WisconsinCalculus/blob/master/figures/221/09surf_of_rotation2.py
# plot surface of revolution around x axis between [0, 3]
# best if r(t) descreases
rad(t) = (t = clamp(t, -1, 1); sqrt(1 - t^2))
rad2(t) = 1/2
viewp = [2,-2,1]
##
unitize(x) = x / norm(x)
"""Orthogonal projection along the vector viewp"""
function make_Pmat(viewp)
a = unitize( [-viewp[2], viewp[1], 0] )
b = unitize( [-viewp[3]*viewp[1],
-viewp[3]*viewp[2],
viewp[1]^2 + viewp[2]^2] )
collect(zip(a,b))
end
#linear projection of R^3 onto R^2
function proj(X, viewp)
Pmat = make_Pmat(viewp)
x=sum([Pmat[i][1]*X[i] for i in 1:3])
y=sum([Pmat[i][2]*X[i] for i in 1:3])
(x, y) # a point
end
proj(X) = proj(X, viewp)
# discrete determinant of Jacobian; area multiplier?
function jac(X, u, v)
ϵ = 0.000001
A = map((p,q) -> (p-q)/ϵ, X(u+ϵ/2, v), X(u-ϵ/2, v))
B = map((p,q) -> (p-q)/ϵ, X(u, v+ϵ/2), X(u, v-ϵ/2))
return A[1]*B[2]-A[2]*B[1]
end
# surface of revolution about the z axis
surf(t, z) = (rad(t)*cos(z), rad(t)*sin(z), t)
surf2(t, z) = (rad2(t)*cos(z), rad2(t)*sin(z), t)
# project the surface at (t, a=theta)
psurf(t,z) = proj(surf(t,z))
psurf2(t, z) = proj(surf2(t,z))
bisect(f, a, b) = find_zero(f, (a,b), Bisection())
_fold(t, zmin, zmax) = bisect(z -> jac(psurf, t, z), zmin, zmax)
_foldz(z, tmin, tmax) = bisect(t -> jac(psurf, t, z), tmin, tmax)
# create shape holding project disc
drawdiscF(t) = Shape([psurf(t, 2*i*pi/100) for i in 1:101])
drawdiscI(t) = Shape([psurf2(t, 2*i*pi/100) for i in 1:101])
# project a line between two points
pline!(p,q; kwargs...) = plot!([proj(p),proj(q)];
line_style..., kwargs...)
α = 1.0
line_style = (; line=(:black, 1))
plot(; empty_style..., aspect_ratio=:equal)
# washer
t0 = sqrt(3/4)
Δ = .03
δ = 0.785398 + 0.05
x₀ = -.25
plot!(drawdiscF(x₀-Δ); fill=(:black,), line=(:black,1))
plot!(drawdiscF(x₀); fill=(:orange,), line=(:black,1))
plot!(drawdiscI(x₀); fill=(:white,1.0), line=(:black,1))
x₀ = 0.35
plot!(drawdiscF(x₀-Δ); fill=(:black,), line=(:black,1))
plot!(drawdiscF(x₀); fill=(:orange,), line=(:black,1))
plot!(drawdiscI(x₀); fill=(:white,1.0), line=(:black,1))
z0 = 3pi/2 - δ
pline!(surf(t0, z0), surf(-t0, z0); line=(:black, 1))
pline!(surf(t0, z0+pi), surf(-t0, z0+pi); line=(:black, 1))
# boundary of sphere
z0 = 3pi/2 - δ
curve = [psurf(t, z0) for t in range(-t0, t0, 100)]
plot!(curve; line=(:black,3))
z0 = 3pi/2 - δ + pi
curve = [psurf(t, z0) for t in range(-t0, t0, 100)]
plot!(curve; line=(:black,3))
# caps
curve = [psurf(t0, θ) for θ in range(0, 2pi, 100)]
plot!(curve, line=(:black, 2))
curve = [psurf(-t0, θ) for θ in range(0, 2pi, 100)]
plot!(curve, line=(:black, 2))
# Shade lines
δ = pi/6
Δₜ = (4pi/2 - (3pi/2 - δ))/(2*25)
for θ ∈ range(3pi/2-δ, 4pi/2, 25)
curve = [psurf(t, θ) for t in
range(-t0, max(-t0, -t0 + 1/2*sin(θ+δ+pi/2 + pi/2)), 20)]
plot!(curve, line=(:black, 1))
curve = [psurf(t, θ+Δₜ) for t in
range(-t0, max(-t0, -t0 + 1/3*sin(θ+δ+pi/2 + pi/2)), 20)]
plot!(curve, line=(:black, 1))
end
#=
f1 = [[t, _fold(t, 0, pi/2)] for t in range(-0.5, -0.1, 26)]
for f in f1
plot!([psurf( f[1], f[2]-k*0.01*(6-f[1]) )
for k in 1:21]; line=(:black, 1))
end
=#
current()
end
plt
```
Figure showing sphere with interior cylinder bored out.
:::
This cross-sectional figure is used to better understand the key dimensions.
::: {#fig-napkin-ring-2}
```{julia}
#| hold: true
#| echo: false
#The following illustrates $R=5$ and $h=8$.
#The following illustrates $R=1$ and $h=2sqrt(3/4$.
plt = let
gr()
empty_style = (xaxis=([], false),
yaxis=([], false),
framestyle=:origin,
legend=false)
R =5; h1 = 2*4
R = 1; h = 2*sqrt(3/4)
theta = asin(h1/2/R)
thetas = range(-theta, stop=theta, length=100)
ts = range(-pi, stop=pi, length=100)
y = h1/4
θ = theta = asin(h/2/R)
thetas = range(-theta, stop=theta, length=100)
ts = range(-pi, stop=pi, length=100)
y = h/4
p = plot(legend=false, aspect_ratio=:equal);
plot!(p, R*cos.(ts), R*sin.(ts));
plot!(p, R*cos.(thetas), R*sin.(thetas), color=:orange);
plot(; empty_style..., aspect_ratio=:equal)
plot!(R*cos.(ts), R*sin.(ts); line=(:black,));
plot!(R*cos.(thetas), R*sin.(thetas), line=(:orange,1));
plot!(p, [R*cos.(theta), R*cos.(theta)], [h1/2, -h1/2], color=:orange);
plot!(p, [R*cos.(theta), sqrt(R^2 - y^2)], [y, y], color=:orange)
plot!([R*cos.(theta), R*cos.(theta)], [h/2, -h/2]; color=:orange);
plot!([R*cos.(theta), sqrt(R^2 - y^2)], [y, y]; line=(:orange,3))
plot!(p, [0, R*cos.(theta)], [0,0], color=:red);
plot!(p,[ 0, R*cos.(theta)], [0,h1/2], color=:red);
plot!([0, R*cos.(theta)], [0,0], color=:red);
plot!([ 0, R*cos.(theta)], [0,h/2], color=:red);
annotate!(p, [(.5, -2/3, "sqrt(R²- (h/2)²)"),
(R*cos.(theta)-.6, h1/4, "h/2"),
(1.5, 1.75*tan.(theta), "R")])
p
x₀ = sqrt(R^2 - (h/2)^2)
annotate!( [
(x₀/2, 0, text(L"\sqrt{R^2- (\frac{h}{2})^2}",10, :top)),
(x₀, h/4, text(L"\frac{h}{2}",:right)),
(R/2*cos(θ),R/2*sin(θ), text(L"R", :bottom; rotation=rad2deg(θ)))
])
current()
end
plt
```
```{julia}
#| echo: false
plotly()
nothing
```
Side view illustrating key dimensions of napkin ring problem with $R$ being the radius of the sphere and $h$ being the height of the resulting interior cylinder.
:::
The small orange line is rotated, so using the washer method we get the cross sections given by $\pi(r_0^2 - r_i^2)$, the outer and inner radii, as a function of $y$.