diff --git a/quarto/adjust_plotly.jl b/quarto/adjust_plotly.jl index ec06770..91bac71 100644 --- a/quarto/adjust_plotly.jl +++ b/quarto/adjust_plotly.jl @@ -9,7 +9,6 @@ function _add_plotly(f) #lineno = 117 - r = readlines(f) inserted = false open(f, "w") do io diff --git a/quarto/integrals/arc_length.qmd b/quarto/integrals/arc_length.qmd index a3d6e5b..e74509a 100644 --- a/quarto/integrals/arc_length.qmd +++ b/quarto/integrals/arc_length.qmd @@ -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 diff --git a/quarto/integrals/area.qmd b/quarto/integrals/area.qmd index 8cbe09d..b5701a1 100644 --- a/quarto/integrals/area.qmd +++ b/quarto/integrals/area.qmd @@ -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$: diff --git a/quarto/integrals/twelve-qs.qmd b/quarto/integrals/twelve-qs.qmd index 4622360..f8a9485 100644 --- a/quarto/integrals/twelve-qs.qmd +++ b/quarto/integrals/twelve-qs.qmd @@ -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() ``` diff --git a/quarto/integrals/volumes_slice.qmd b/quarto/integrals/volumes_slice.qmd index 4b4d0d7..80d466f 100644 --- a/quarto/integrals/volumes_slice.qmd +++ b/quarto/integrals/volumes_slice.qmd @@ -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$.