CalculusWithJulia.jl

CalculusWithJulia

A package to accompany notes at https://calculuswithjulia.github.io on using Julia for topics from the calculus sequence.

This package does two things:

• It loads a few other packages making it easier to use (and install) the functionality provided by them and

• It defines a handful of functions for convenience. The exported ones

are unzip, rangeclamp tangent, secant, D (and the prime notation), divergence, gradient, curl, and ∇, along with some plotting functions. The constant e is assigned to exp(1).

Packages loaded by CalculusWithJulia

• The SpecialFunctions is loaded giving access to a few special functions used in these notes, e.g., airyai, gamma

• The ForwardDiff package is loaded giving access to its derivative, gradient, jacobian, and hessian functions for finding automatic derivatives of functions. In addition, this package defines ' (for functions) to return a derivative (which commits type piracy), ∇ to find the gradient (∇(f)), the divergence (∇⋅F). and the curl (∇×F), along with divergence and curl.

• The LinearAlgebra package is loaded for access to several of its functions for working with vectors norm, cdot (⋅), cross (×), det.

• The PlotUtils package is loaded so that its adapted_grid function is available.

Packages with extra features added when loaded

Either through extensions or through the Julia package Requires there is additional code to be run when the following packages load:

• SymPy: for symbolic math.

• Plots: the Plots package provides a plotting interface.

Several plot recipes are provided to ease the creation of plots in the notes. plotif, trimplot, and signchart are used for plotting univariate functions; plot_polar and plot_parametric are used to plot curves in 2 or 3 dimensions; plot_parametric also makes the plotting og parameterically defined surfaces easier; vectorfieldplot and vectorfieldplot3d can be used to plot vector fields; and arrow is a simplified interface to quiver that also indicates 3D vectors.

The plot_implicit function can plot 2D implicit plots. (It is borrowed from ImplicitPlots.jl, which is avoided, as it has dependencies that hold other packages back.)

Other packages with a recurring role in the accompanying notes:

• Roots is used to find zeros of univariate functions

• SymPy for symbolic math

• QuadGK and HCubature are used for numeric integration

SignChart(f, a, b)

Numerically identifies values of x in [a,b] where f is 0, oo, or undefined.

Displays the output as a sideways interval displaying the sign of f in between these values.

Example

julia> SignChart((x -> sqrt(1 - x^2))', -1, 1)
         ↑
         ⋮
        1.0         is infinite
         ⋮
         +
         ⋮
        0.0         a zero
         ⋮
         -
         ⋮
       -1.0         is infinite
         ⋮
         ↓

D(f)

Function interface to ForwardDiff.derivative.

A method for adjoint for functions dispatches to D, so that the notation f' can be used to take the derivative of a function. (This is type piracy.)

arrow!(p, v)

Add the vector v to the plot anchored at p.

This would just be a call to quiver, but there is no 3-D version of that. As well, the syntax for quiver is a bit awkward for plotting just a single arrow. (Though efficient if plotting many).

using Plots
r(t) = [sin(t), cos(t), t]
rp(t) = [cos(t), -sin(t), 1]
plot(unzip(r, 0, 2pi)...)
t0 = 1
arrow!(r(t0), rp(t0))

curl(F)

Find curl of a 2 or 3-D vector field.

divergence(F)

Find divergence of a 3-D vector vield.

fisheye(f)

Transform f defined on (-∞, ∞) to a new function whose domain is in (-π/2, π/2) and range is within (-π/2, π/2). Useful for finding all zeros over the real line. For example

f(x) = 1 + 100x^2 - x^3
fzeros(f, -100, 100) # empty just misses the zero found with:
fzeros(fisheye(f), -pi/2, pi/2) .|> tan  # finds 100.19469143521222, not perfect but easy to get

By Gunter Fuchs.

fubini(f, [zs], [ys], xs; rtol=missing, kws...)

Integrate f of 1, 2, or 3 input variables.

The zs may depend (x,y), the ys may depend on x

Examples

# integrate over the unit square
fubini((x,y) -> sin(x-y), (0,1), (0,1))

# integrate over a triangle
fubini((x,y) -> 1, (0,identity), (0,1 ))

#
f(x,y,z) = x*y^2*z^3
fubini(f, (0,(x,y) ->  x+ y), (0, x -> x), (0,1))

!!! Note This uses nested calls to quadgk. The use of hcubature is recommended, typically after a change of variables to make a rectangular domain. The relative tolerance increases at each nested level.

lim(f, c; n=6, m=1, dir="+-")
lim(f, c, dir; n-5)

Means to generate numeric table of values of f as h gets close to c.

• n, m: powers of 10 to add (subtract) to (from) c.
• dir: Either "+-" (show left and right), "+" (right limit), or "-" (left limit). Can also use functions +, -, ±.

Example:

julia> f(x) = sin(x) / x
f (generic function with 1 method)

julia> lim(f, 0)
 0.1        0.9983341664682815
 0.01       0.9999833334166665
 0.001      0.9999998333333416
 0.0001     0.9999999983333334
 1.0e-5     0.9999999999833332
 1.0e-6     0.9999999999998334
   ⋮          ⋮
   c          L?
   ⋮          ⋮
-1.0e-6     0.9999999999998334
-1.0e-5     0.9999999999833332
-0.0001     0.9999999983333334
-0.001      0.9999998333333416
-0.01       0.9999833334166665
-0.1        0.9983341664682815

newton_plot!(f, x0; steps=5, annotate_steps::Int=0, kwargs...)

Add trace of Newton's method to plot.

• steps: how many steps from x0 to illustrate
• annotate_steps::Int: how may steps to annotate

newton_vis(f, x0, a=Inf, b=-Inf; steps=5, kwargs...)

Visualize `F(x,y,z) = c` by plotting assorted contour lines

This graphic makes slices in the x, y, and/or z direction of the 3-D level surface and plots them accordingly. Which slices (and their colors) are specified through a dictionary.

Examples:

F(x,y,z) = x^2 + y^2 + x^2
plot_implicit_surface(F, 20)  # 20 slices in z direction
plot_implicit_surface(F, 20, slices=Dict(:x=>:blue, :y=>:red, :z=>:green), nlevels=6) # all 3 shown

# A heart
a,b = 1,3
F(x,y,z) = (x^2+((1+b)*y)^2+z^2-1)^3-x^2*z^3-a*y^2*z^3
plot_implicit_surface(F, xlims=-2..2,ylims=-1..1,zlims=-1..2)

Note: Idea from.

Not exported.

plot_parametric(ab, r; kwargs...)
plot_parametric!(ab, r; kwargs...)
plot_parametric(u, v, F; kwargs...)
plot_parametric!(u, v, F; kwargs...)

Make a parametric plot of a space curve or parametrized surface

The intervals to plot over are specifed using a..b notation, from IntervalSets

plotif(f, g, a, b)

Plot of f over [a,b] with the intervals where g ≥ 0 highlighted in many ways.

rangeclamp(f, hi=20, lo=-hi; replacement=NaN)

Modify f so that values of f(x) outside of [lo,hi] are replaced by replacement.

Examples

f(x) = 1/x
plot(rangeclamp(f), -1, 1)
plot(rangeclamp(f, 10), -1, 1) # no `abs(y)` values exceeding 10

riemann(f, a, b, n; method="right"

Compute an approximations to the definite integral of f over [a,b] using an equal-sized partition of size n+1.

method: "right" (default), "left", "trapezoid", "simpsons", "ct", "m̃" (minimum over interval), "M̃" (maximum over interval)

Example:

f(x) = exp(x^2)
riemann(f, 0, 1, 1000)   # default right-Riemann sums
riemann(f, 0, 1, 1000; method="left")       # left sums
riemann(f, 0, 1, 1000; method="trapezoid")  # use trapezoid rule
riemann(f, 0, 1, 1000; method="simpsons")   # use Simpson's rule

riemann_plot!(f, a, b, n; method="method", fill, kwargs...)
riemann_plot(f, a, b, n; method="method", fill, kwargs...)

Add visualization of riemann sum in a layer.

• method: one of right, left, trapezoid, simpsons
• fill: to specify fill color, something like ("green", 0.25, 0) will fill in green with an alpha transparency.

secant(f::Function, a, b)

Returns a function describing the secant line to the graph of f at x=a and x=b.

Example. Where does the secant line intersect the y axis?

f(x) = sin(x)
a, b = pi/4, pi/3
sl = secant(f, a, b)  # or sl(x) = secant(f, a, b)(x) to use a generic function
sl(0)

sign_chart(f, a, b; atol=1e-4)

Create a sign chart for f over (a,b). Returns a collection of named tuples, each with an identified zero or vertical asymptote and the corresponding sign change. The tolerance is used to disambiguate numerically found values.

Example

julia> sign_chart(x -> (x-1/2)/(x*(1-x)), 0, 1)
3-element Vector{NamedTuple{(:zero_oo_NaN, :sign_change)}}:
 (zero_oo_NaN = 0.0, sign_change = an endpoint)
 (zero_oo_NaN = 0.5, sign_change = - to +)
 (zero_oo_NaN = 1.0, sign_change = an endpoint)

signchart(f, a, b)

Plot f over a,b with different color when negative.

tangent(f::Function, c)

Returns a function describing the tangent line to the graph of f at x=c.

Example. Where does the tangent line intersect the y axis?

f(x) = sin(x)
tl = tangent(f, pi/4)  # or tl(x) = tangent(f, pi/3)(x) to use a generic function
tl(0)

Uses the automatic derivative of f to find the slope of the tangent line at x=c.

trimplot(f, a, b, c=20; kwargs...)

Plot f over [a,b] but break graph if it exceeds c in absolute value.

unzip(vs)
unzip(v1, v2, ...)
unzip(r::Function, a, b)

Take a vector of points described by vectors (as returned by, say r(t)=[sin(t),cos(t)], r.([1,2,3]), and return a tuple of collected x values, y values, and optionally z values.

Wrapper around the invert function of SplitApplyCombine.

If the argument is specified as a comma separated collection of vectors, then these are combined and passed along.

If the argument is a function and two end points, then the function is evaluated at points between a and b; for the univaraite case, the points are chosen adaptively.

This is useful for plotting when the data is more conveniently represented in terms of vectors, but the plotting interface requires the x and y values collected.

Examples:

using Plots
r(t) = [sin(t), cos(t)]
rp(t) = [cos(t), -sin(t)]
plot(unzip(r, 0, 2pi)...)  # calls plot(xs, ys)

t0, t1 = pi/6, pi/4

p, v = r(t0), rp(t0)
plot!(unzip(p, p+v)...)  # connect p to p+v with line

p, v = r(t1), rp(t1)
quiver!(unzip([p])..., quiver=unzip([v]))

Based on unzip from the Plots package. Implemented through invert of SplitApplyCombine

Note: for a vector of points, xs, each of length 2, a similar functionality would be (first.(xs), last.(xs)). If each point had length 3, then with second(x)=x[2], a similar functionality would be (first.(xs), second.(xs), last.(xs)).

```

uvec(x)

Helper to find a unit vector.

vectorfieldplot(V; xlim=(-5,5), ylim=(-5,5), n=10; kwargs...)

V is a function that takes a point and returns a vector (2D dimensions), such as V(x) = x[1]^2 + x[2]^2.

The grid xlim × ylim is paritioned into (n+1) × (n+1) points. At each point, pt, a vector proportional to V(pt) is drawn.

This is written to add to an existing plot.

plot()  # make a plot
V(x,y) = [x, y-x]
vectorfield_plot!(p, V)
p