CalculusWithJuliaSquared API

CalculusWithJuliaSquared

A personal, pure-Julia fork of CalculusWithJulia.jl 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).

• It supplies the exact symbolic algebra that Symbolics core does not: exact_trig_values

(the special-angle table, so cos(π/6) becomes √3/2 instead of staying unevaluated), factored_poly and poly_factors (factoring over the rationals), and partial_fractions. These stand in for SymPy's automatic special angles, factor and apart. Alongside them, numeric_roots and root_enclosures find roots exactly and report them as floats or as certified intervals, standing in for N.(solve(...)) and sympy.real_roots; and divrem divides one polynomial by another.

Packages loaded by CalculusWithJuliaSquared

• 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.

• The Symbolics package is loaded (and reexported) giving access to symbolic math (@variables, etc.) along with symbolic gradient, divergence, and curl methods – pure Julia, no Python dependency.

• The Nemo package is loaded – imported, not reexported – which switches on Symbolics.symbolic_solve for polynomial equations, and also backs factored_poly, poly_factors and partial_fractions. The module name itself IS exported, so Nemo.overlaps, Nemo.midpoint and Nemo.radius are available for the balls root_enclosures returns; but no using Nemo is needed downstream, and none of Nemo's own names (derivative, coeff, roots, ...) enter the namespace, where they would collide with Symbolics.

• The Plots package is loaded (and reexported) providing the plotting interface directly – no separate using Plots needed.

• The LaTeXStrings package is loaded (and reexported), so L"..." works with no separate using. Plots does not pass this through. It is carried here because it is house standard across the sibling study repos, which means a downstream package or notebook environment need not name it alongside this one.

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

• QuadGK and HCubature are used for numeric integration

Cautions for anyone else using this package

This is a personal study fork, not a registered package: you have to go out of your way to install it, and nobody maintains it for you. Several of its conveniences change Julia's behaviour globally for the whole session, not only for calls into this package. They are collected here so that nobody has to discover them by debugging.

Type piracy: what this package does, and the rule it follows

Type piracy means adding a method to a function you do not own, dispatching on a type you do not own. It matters because Julia's method tables are global: such a method is visible to every package in the session, not only to code that calls into this one. Loading this package can therefore change how unrelated code behaves, which is why the practice is discouraged.

This package commits it five times, each deliberately, and each passing the same test:

The benign test. Every call the pirated method answers is a call that would otherwise have thrown — a MethodError, or an outright error. No working code changes its behaviour; code that used to fail now succeeds. The manual's own carve-out for tightly coupled packages that "separate features from definitions" is the ground these stand on.

A piracy that changed the result of a call that already worked would not pass that test, and none of the five below does.

Three of these are worth a further word:

• divrem is named, not renamed. A poly_divrem of our own would have avoided the piracy entirely. It is not used because the point being taught is that Julia's generic divrem divides polynomials exactly as it divides integers; a bespoke name would state the opposite.

• Plots._show pirates an internal. The leading underscore marks it as private to Plots, so unlike the other four it carries no API stability promise at all: a patch release could rename or remove it, and the symptom would be a plot that silently stops being interactive rather than an error. Measured 2026-09-11 against Plots v1 — _best_html_output_type maps :plotly => :html, and the generic _show(::IO, ::MIME"text/html", ::Plot) has no :html branch, so without this method the call throws. That is what keeps it on the benign side of the test above.

• find_zero may one day collide. Should Roots add its own Symbolics support, expect a method-overwrite warning on load. That is not a bug to work around: the fix is to delete our block, because upstream's version supersedes it. The same applies to any of the five if the owning package adopts the method itself.

Nemo is imported, and only its name is exported

Nemo does two jobs here: it switches on Symbolics.symbolic_solve for polynomial equations (via Symbolics' SymbolicsNemoExt), and it backs factored_poly, poly_factors, partial_fractions, numeric_roots, root_enclosures and divrem. So loading this package changes what Symbolics itself can do — code that fails without it will succeed with it.

How it is loaded is a deliberate middle course, and knowing which one you are in explains every "why must I qualify this?" question:

The module name is exported, exactly as ForwardDiff's is, so you can call Nemo.overlaps, Nemo.midpoint, Nemo.radius and Nemo.contains on the balls root_enclosures returns without adding Nemo to your own project. No using Nemo is needed, and none of Nemo's functions enter your namespace.

Why not reexport. Nemo exports roots, degree, derivative, coeff, factor, term and terms. Every one of those collides with something these notes use — Polynomials.roots and Polynomials.degree above all, and derivative is a three-way clash between Nemo, Polynomials and Symbolics. Reexporting would turn working code into ambiguity errors. The same reasoning applies to Symbolics' own public-but-unexported names (derivative, value, get_variables, jacobian, hessian and others): this package does not re-export them either, so they stay qualified as Symbolics.derivative. That is upstream's decision and overriding it would break the collision-avoidance it exists for.

A lot of names arrive at once

Roots, LinearAlgebra, SpecialFunctions, IntervalSets, Symbolics, Plots and LaTeXStrings are reexported; ForwardDiff and Nemo are exported as module names; and e is exported as exp(1). Clashes are real rather than theoretical: alongside SciML's BracketingNonlinearSolve, both Bisection and solve become ambiguous and have to be qualified.

The plotting backend is configured on load. __init__ selects GR and forces headless mode whenever Julia is non-interactive, so that document renders embed figures instead of trying to open a window.

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
         ⋮
         ↓

divrem(a::Num, b::Num)

Divide one polynomial by another, returning (quotient, remainder).

This is Julia's generic divrem extended to symbolic polynomials, and it means what it means for integers: a == b*q + r, with the remainder of strictly lower degree than b. It is the division algorithm behind rewriting a rational expression as a polynomial plus a proper fraction, which is how a slant asymptote is found.

julia> using CalculusWithJuliaSquared

julia> @variables x;

julia> q, r = divrem(5x^3 + 6x^2 + 2, x - 1)
(11 + 11x + 5(x^2), 13)

Both arguments must be polynomials over the rationals in the same single variable, which is read off the expressions themselves since there is no argument naming it. A constant on one side is fine. Dividing by zero throws DivideError, as it does for numbers.

For the partial-fraction decomposition of the whole rational expression in one step, see partial_fractions.

show(io, ::MIME"text/latex", x::Symbolics.Num)
show(io, ::MIME"text/html",  x::Symbolics.Num)

Typeset a symbolic expression as display math in HTML/LaTeX frontends — Quarto, Jupyter, Documenter — parallel to the text/latex show SymPy has built in. Companion methods do the same for the vector of roots Symbolics.symbolic_solve returns, which would otherwise print its full internal type name ahead of the mathematics.

Latexify's own text/latex output wraps the expression in $$\begin{equation}...\end{equation}$$, which Quarto renders literally, so these emit a clean \[ ... \] through both MIME types instead. There is no effect in the plain-text REPL, which uses text/plain.

D(f)
D(f, n::Int)

Function interface to ForwardDiff.derivative; D(f, n) applies it n times, with D(f, 0) returning f itself.

A method for Base.adjoint on functions dispatches to D, so that the notation f' can be used to take the derivative of a function.

julia> using CalculusWithJuliaSquared

julia> f(u) = u^2;

julia> f'(3.0)
6.0

What the samples on ONE side of c say. s is +1 for the right, -1 for the left.

• nothing — ex cannot be evaluated on that side, so its domain does not reach c from that direction. Not a failure: log at 0 is exactly this.
• (:finite, y) — the samples are settling on y.
• (:diverges, ±Inf) — the increments are not collapsing and hold one sign.
• (:erratic, NaN) — evaluable, but neither settling nor monotone. No opinion.

Testing the increments rather than the magnitude is what catches logarithmic divergence, which grows by a constant per decade and never passes a fixed threshold.

Is a substituted result usable as a limit?

A numeric result must be finite. A symbolic result is fine too — 5x^4 is the honest answer to a limit taken in h — provided it no longer mentions the limit variable and grounds to a finite value. Requiring a Number here was a bug: it threw away correct answers for every limit carrying a parameter, and sent them to the Gruntz engine instead.

arrow(p, v)

Add vector, v, to plot anchored at point p.

Example

Fn = parametric(t -> [2cos(t), 3sin(t)])
Fnp = t -> ForwardDiff.derivative(Fn, t)
p = plot(Fn, 0, 2pi, legend=false)
for t in 0:pi/4:pi
   arrow!(Fn(t), Fnp(t))
end
p

arrow(p, v)

Add vector, v, to plot anchored at point p.

Example

Fn = parametric(t -> [2cos(t), 3sin(t)])
Fnp = t -> ForwardDiff.derivative(Fn, t)
p = plot(Fn, 0, 2pi, legend=false)
for t in 0:pi/4:pi
   arrow!(Fn(t), Fnp(t))
end
p

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

arrow(p, v)

Add vector, v, to plot anchored at point p.

Example

Fn = parametric(t -> [2cos(t), 3sin(t)])
Fnp = t -> ForwardDiff.derivative(Fn, t)
p = plot(Fn, 0, 2pi, legend=false)
for t in 0:pi/4:pi
   arrow!(Fn(t), Fnp(t))
end
p

conventional_latex(ex) -> String

Typeset a symbolic expression the way a mathematics text writes it.

Latexify renders a Symbolics expression faithfully, and faithfully means printing the canonical algebraic form rather than conventional notation. This walks the expression and emits the LaTeX directly, fixing three things at once:

It also removes a genuine defect rather than an infelicity: a leading negative renders with a space after the opening delimiter, and Pandoc will not open inline math on a dollar-sign-then-space, so those cells appear on the page as a literal dollar followed by raw LaTeX.

The result is the body of a math expression, with no delimiters, so a caller wraps it as it needs – inline, display, or as a cell of a table.

julia> @variables x;

julia> conventional_latex(x - 1)
"x - 1"

julia> conventional_latex((2//3) * Symbolics.Num(pi))
"\\frac{2 \\pi}{3}"

Sums are ordered by descending total degree, so the constant comes last, except that a term which is nothing but a radical sorts last of all – which is how -b + \sqrt{b^2-4c} and 3 + \sqrt{41} are conventionally written. The ordering is total, so a rebuild renders a given expression identically; it will not quietly rearrange a published page.

When every term of a sum is a fraction over the same denominator they are written over one bar, which is what turns the pieces above into the quadratic formula rather than two fractions added together. Terms over different denominators are left alone – otherwise a partial-fraction decomposition, whose entire point is separate denominators, would be recombined into the thing it was decomposed from.

Any expression shape this does not recognise is passed to Latexify unchanged, so an unfamiliar function renders as it always did rather than failing.

curl(F)

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

divergence(F)

Find divergence of a 3-D vector vield.

exact_trig_values(ex)

Replace every trigonometric function applied to a rational multiple of π in ex with its exact value, leaving the rest of the expression alone.

Symbolics has no table of special angles, so cos(Num(π)/6) simply stays unevaluated, and simplify makes matters worse by folding the π to a float. This walks the expression instead and substitutes the exact value, which is an ordinary symbolic term: sqrt(3)/2 prints as √3/2 and typesets as such.

Handles sin, cos, tan, csc, sec and cot at any rational multiple of π with denominator 1, 2, 3, 4 or 6, in every quadrant.

Anything else is returned unchanged – an angle that is not such a multiple, an argument still carrying a free variable, and the poles (tan(π/2), cot(0)), which have no value to give. Note that the angle has to be exact going in: a Float64 that merely rounds to π/6 is refused rather than guessed at.

Examples

julia> PI = Symbolics.Num(pi);

julia> exact_trig_values(cos(PI/6))
sqrt(3) / 2

julia> exact_trig_values.(cos.([0, PI/6, PI/4, PI/3, PI/2]))
5-element Vector{Num}:
         1
 sqrt(3) / 2
 sqrt(2) / 2
       1//2
         0

julia> exact_trig_values(cos(Symbolics.Num(0.5235987755982988)))   # a float, not π/6
cos(0.5235987755982988)

See also factored_poly, partial_fractions.

factored_poly(ex, var)

Factor the univariate polynomial ex over the rationals and return the factored expression.

This is poly_factors multiplied back together, and carries the same restrictions: the factorisation is over the rationals, so x^2 - 2 is returned unchanged, and a non-polynomial or non-rational coefficient throws.

Examples

julia> @variables x;

julia> factored_poly(x^3 - 6x^2 + 11x - 6, x)
(-3 + x)*(-2 + x)*(-1 + x)

julia> factored_poly(x^2 - 2, x)
-2 + x^2

See also poly_factors, partial_fractions.

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
find_zeros(f, -100, 100) # empty just misses the zero found with:
find_zeros(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

numeric_roots(ex, var; real_only=false)

Every root of the univariate polynomial ex, as a floating-point number.

Returns ComplexF64 values by default, or Float64 values when real_only=true, in which case the non-real roots are dropped. This is the counterpart of SymPy's N.(solve(...)) and, with real_only=true, of sympy.real_roots.

Roots are found exactly, as algebraic numbers, and rounded only on the way out. So they are returned even for polynomials that have no formula in radicals, where symbolic_solve can only answer roots_of(...):

julia> using CalculusWithJuliaSquared

julia> @variables x;

julia> numeric_roots(x^2 - 2, x; real_only=true)
2-element Vector{Float64}:
 -1.4142135623730951
  1.4142135623730951

Two guarantees worth relying on when rendering a document:

• The order is stable. Real results are sorted ascending; complex results are sorted by real part, then imaginary part. The same input renders the same way every time.
• Repeated roots repeat. A double root appears twice, so with real_only=false the number of values always equals the degree.

ex must be a polynomial in var alone with rational coefficients; anything else throws rather than guessing. A constant is refused, since it has either no roots or all of them.

For real roots with an error bound you can reason about – rather than a float whose accuracy is unstated – use root_enclosures.

See also poly_factors, factored_poly.

partial_fractions(ex, var)

Decompose the rational expression ex into partial fractions over the rationals, returning the decomposition as a symbolic sum (SymPy spells this apart).

ex must be a single quotient of polynomials in var – the shape p/q – or a polynomial, which is returned unchanged. A sum of separate fractions is not recognised; put it over a common denominator first with simplify or simplify_fractions.

The polynomial part of an improper fraction is included in the sum, and repeated factors in the denominator produce the expected higher-power terms.

Examples

julia> @variables x;

julia> partial_fractions(1/((x-1)*(x-2)), x)
-1 / (-1 + x) + 1 / (-2 + x)

julia> partial_fractions((x+3)/((x-1)^2*(x+2)), x)
(-1//9) / (-1 + x) + (1//9) / (2 + x) + (4//3) / ((-1 + x)^2)

julia> partial_fractions(x^3/((x-1)*(x-2)), x)      # improper: polynomial part included
3 + x + -1 / (-1 + x) + 8 / (-2 + x)

See also factored_poly, exact_trig_values.

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.

poly_factors(ex, var)

Factor the univariate polynomial ex over the rationals and return its factors as a vector, each repeated according to its multiplicity, with any constant factor first.

Use this where a count is wanted – length(poly_factors(ex, x)) – and factored_poly where the factored expression itself is wanted.

Factoring is over the rationals, so x^2 - 2 comes back as a single factor: the factorisation (x-√2)(x+√2) exists but is not rational. symbolic_solve will give those roots.

The order is sorted by degree, so repeated runs and repeated renders agree.

Throws if ex is not a polynomial in var alone, or if any coefficient is not rational.

Examples

julia> @variables x;

julia> poly_factors(x^3 - 6x^2 + 11x - 6, x)
3-element Vector{Num}:
 -3 + x
 -2 + x
 -1 + x

julia> length(poly_factors(x^12 - 1, x))
6

julia> poly_factors(x^2 - 2, x)          # irreducible over the rationals
1-element Vector{Num}:
 -2 + x^2

See also factored_poly, partial_fractions.

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.

root_enclosures(ex, var; bits=128)

The real roots of the univariate polynomial ex, each as an interval that is guaranteed to contain one.

Returns Nemo.ArbFieldElem values – ball arithmetic, printed as [midpoint +/- radius] – sorted ascending. Each is a rigorous enclosure computed from the exact algebraic root, not a float with an informal error estimate.

The point of an enclosure is that it supports proof rather than eyeballing. Two nearby floats cannot tell you whether there are two roots or one root computed twice; two enclosures that do not overlap can only come from two distinct roots:

julia> @variables s;

julia> es = root_enclosures(s^15 - 16129s^2 + 254s - 1, s; bits=96);

julia> Nemo.overlaps(es[1], es[2])      # disjoint => genuinely two roots
false

bits sets the working precision, and is a real dial rather than decoration. On the polynomial above the first two enclosures still overlap at 53 bits – the honest answer being "these cannot be told apart yet" – and separate only at 64. That threshold is the point of the whole exercise: 53 bits is exactly the precision of a Float64 mantissa, so no computation carried in double precision can establish that this cluster is two roots rather than one. Nothing built on Float64 interval endpoints could do it either, since their radius bottoms out at an ulp.

Useful Nemo functions for working with the results – Nemo is imported by this package but not reexported, so they must be qualified:

• Nemo.overlaps(a, b) – do two enclosures intersect?
• Nemo.contains(ball, x) – is x inside? Note: there is no method for a Float64 argument; wrap it first, as Nemo.contains(b, Nemo.ArbField(96)(6.94)).
• Nemo.midpoint(ball), Nemo.radius(ball) – the two halves of [m +/- r].

Arithmetic works and stays certified (b + 1.0, b^2, sin(b), sqrt(b)), and Float64.(root_enclosures(...)) drops back to plain numbers when that is all you need.

Only real roots are returned; a ball is real by construction. For the complex ones, or for plain floats, use numeric_roots.

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.

symlim(ex, v, c; side = :both, cancel = true, check = true, n = 8, secs = 10)

Symbolic limit of the expression ex as the variable v approaches c.

Returns (value, route), where route names the method that produced the answer. A value of nothing means every method declined — an honest refusal rather than a guess.

Routes, in the order they are tried

julia> @variables x::Real;

julia> symlim((x^2 - 1)/(x - 1), x, 1)
(2, :cancel)

julia> symlim(sin(x)/x, x, 0)
(1//1, :series)

julia> symlim(log(x)/x, x, Inf)
(0, :gruntz)

julia> symlim(x * sin(1/x), x, 0)
(0.0, :squeeze)

julia> symlim(abs(x)/x, x, 0)
(nothing, :sides_disagree)

julia> symlim(abs(x)/x, x, 0; side = :right)
(1//1, :series)

julia> symlim(abs(x)/x, x, 0; side = :left)
(-1//1, :series)

julia> symlim(floor(x), x, 0; side = :right), symlim(floor(x), x, 0; side = :left)
((0, :substitution), (-1, :squeeze))

julia> @variables c::Real;

julia> symlim(3x^2 + c, x, 0; side = :right)
(c, :substitution)

x·sin(1/x) has the limit 0 by the squeeze theorem, which the :squeeze route establishes: interval arithmetic bounds the function on a shrinking neighbourhood of 0, and the enclosures collapse to a point. Where they do not collapse the route declines — sin(x) at infinity encloses to [-1, 1] at every scale, correctly, because it has no limit. The boxes exclude c itself, since a limit never consults f(c); that is what lets floor resolve from the left, where the value at 0 is not the limit.

:substitution is the definition of continuity, computed: the value at c exists and the function approaches it. The route says so. floor from the right is :substitution and from the left is not, which is the statement that floor is right-continuous at the integers.

A free parameter rides through the exact routes — c, 5x^4, 1/x are all honest answers to limits taken in another variable. A route that can only produce a number (:squeeze, :divergent_numeric, the engine at a finite point) is not allowed to invent a value for one, and refuses with :parameter_dependent instead.

Sidedness

A limit exists at c only if the left and right limits both exist and agree, so side = :both refuses when they demonstrably do not: abs(x)/x and 1/x at 0 are jumps, not limits. Ask for :left or :right to get the one-sided answer.

Refusing needs positive evidence from both sides. Where ex simply does not reach c from one direction — log(x) at 0, x^x at 0, anything at the edge of its domain — the defined side is the answer, which is the ordinary reading of lim_{x→0} log(x) = -∞, and every route is asked for that side's limit. Only a genuine two-sided disagreement is refused.

Why the ordering matters

Symbolics computes no limits itself; the available engine, SymbolicLimits, implements the Gruntz algorithm for log-exponential asymptotics at infinity. It is excellent at that and unreliable elsewhere. As of v1.1.5 it does not cancel common factors, so limit((x^2-1)/(x-1), x, 1) returns 0 rather than 2 — with its assumption set reporting confidence. Cancelling and series comparison are tried first because they are both exact and dependable; the engine is asked last, and only for the work it was built for.

Keyword arguments

• cancel = true — clear removable singularities before the engine ever sees the expression. Setting false skips that stage only; the series stage below it will usually still find the right answer, so this is not a way to observe the underlying engine misbehaving. Use check = false for that, or call SymbolicLimits.limit directly.
• side = :both — :left or :right for a one-sided limit. See Sidedness.
• check = true — require every route's answer to agree with the numeric evidence on the side being approached, and move to the next route when it does not. Leave this on: the failure modes it guards against are silent, so nothing else will catch them. Setting false returns the first route's answer unchecked, which is the way to watch the engine misbehave.
• n — highest order used by the series route.
• secs — deadline for the Gruntz stage, which is the only one that can fail to terminate. Start Julia with -t 2 or more, or the watchdog cannot run: a hung call never yields, so a same-thread timer would never fire.

Known limits

• Expansion points involving π. Symbolics.taylor converts π to a float and then to a rational, so a series about π/2 carries a spurious constant term near 1e-17 instead of an exact 0, which defeats leading-order ranking. A Num limit point such as Num(pi)/2 is folded to a float at the door, so it behaves exactly like pi/2 and does not help either. Use lim for those.
• (1 + 1/x)^x as x → ∞. SymbolicLimits does not terminate on this, nor on the log/exp rewrite its own error message suggests; the deadline returns :unresolved.
• Oscillation without a limit. sin(x) as x → ∞ genuinely has none. The :squeeze route treats its enclosure [-1, 1] as a refusal rather than an answer; call IntervalArithmetic directly to see the bound itself, which is the closest analogue to what a SymPy user gets from AccumBounds.
• :squeeze is tried last, because interval arithmetic cannot see that two occurrences of x are the same number: sin(x)/x over [h/10, h] encloses to roughly [0.1, 10] however small h is, and x·floor(1/x) declines the same way. Every exact route is asked first for that reason. The enclosures must close at the rate the scales do, so sqrt(x)·sin(1/x) — whose limit is 0, approached too slowly — declines. Its answer is typed by what the route knows: an enclosure of zero width is exact and reports an integer (floor from the left is -1); one containing zero reports 0.0, not -1.0e-14; anything else is the float midpoint. A float from this route means bounded, not derived — (2/3)(1 - (1/4)^(n+1)) at infinity is 0.6666… here and 2//3 from the Gruntz engine after the a^m → e^{m·log a} rewrite, and the difference is the point.
• log divergence cannot be delegated. SymbolicLimits.limit(log(u), u, 0, :right) returns 0 rather than -Inf (v1.1.5). That answer now fails the check comparison and is discarded, and the numeric increment test supplies the -Inf — load-bearing, not a convenience.
• Unknown symbolic exponents. The order of x^k is k, so the limit of sin(sin(x^2))/x^k at 0 is 0, 1 or ∞ depending on k alone. The series route declines rather than picking one; substitute a concrete k, or rewrite a^x as exp(x·log(a)) where the exponent is no longer the unknown.
• Limits at infinity are mostly the Gruntz engine. Ratios of polynomials go through :reciprocal and come back exact. Everything else is a time-boxed call into the engine, with no numeric evidence to check it against and no series to take. It handles pure power/log/exp forms; add a sin or a sqrt — exp(-x)·sin(x), x/sqrt(x^2+4) — and the answer is :unresolved.

See also tlim, lim.

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.

tlim(num, den, v, c = 0; n = 6)

Limit of num/den as v → c, evaluated by Taylor series.

Both parts are expanded about c and compared by leading order: when the two series start at the same power the limit is the ratio of their leading coefficients; when the numerator starts higher the limit is 0. The answer comes back as an exact Rational wherever the coefficients are exact.

This reaches the limits SymbolicLimits declines outright — anything involving trigonometric functions or roots — because a series turns them into polynomials:

julia> @variables x::Real;

julia> tlim(sin(x), x, x)
1//1

julia> tlim(1 - cos(x), x^2, x)
1//2

julia> tlim(2sin(x) - sin(2x), x - sin(x), x)
6//1

Returns nothing when the method does not apply: an essential singularity, a coefficient that stays symbolic, or a numerator vanishing slower than the denominator (a pole rather than a limit).

Pass side = :left or :right where the expression contains an abs. A series cannot see a sign change: Symbolics.taylor(abs(w), w, 0:n) returns w, which is the right answer only from the right. Given a side, each abs is first resolved against the sign its argument actually holds there, so abs(x)/x at 0 expands to x/x on the right and -x/x on the left:

julia> @variables x::Real;

julia> tlim(abs(x), x, x, 0; side = :right), tlim(abs(x), x, x, 0; side = :left)
(1, -1)

Where the sign is not settled the abs stays put and the method declines.

Two cautions. The expansion point must be one the series can be taken about, so v → ∞ is out of reach — use symlim, which delegates those to the Gruntz engine. And a series argument is circular if used to derive the derivatives its own coefficients assume: proving [sin(x)]' = cos(x) from the Taylor series of sin assumes the answer. Computing lim sin(x)/x is not circular in that way, since the limit is the goal rather than a step toward it.

See also symlim, lim.

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(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=8], [ny=8])

Create a vector field plot using a grid described by xlim, ylim with nx and ny grid points in each direction.

F(x,y) = [-y, x]
vectorfieldplot(F, xlim=(-4,4), ylim=(-4,4))

vectorfieldplot(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=8], [ny=8])

Create a vector field plot using a grid described by xlim, ylim with nx and ny grid points in each direction.

F(x,y) = [-y, x]
vectorfieldplot(F, xlim=(-4,4), ylim=(-4,4))

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

vectorfieldplot(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=8], [ny=8])

Create a vector field plot using a grid described by xlim, ylim with nx and ny grid points in each direction.

F(x,y) = [-y, x]
vectorfieldplot(F, xlim=(-4,4), ylim=(-4,4))

vectorfieldplot3d(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=5], [ny=5])

Create a 3 dimensional vector field plot using a grid described by xlim, ylim, zlim with nx, ny, and nz grid points in each direction.

Note: the vectors are represented with line, not arrow due to no implementation of :quiver3d.

F(x,y,z) = [-y, x,z]
vectorfieldplot3d(F, xlims=(-4,4), ylims=(-4,4), zlims=(0,3))

vectorfieldplot3d(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=5], [ny=5])

Create a 3 dimensional vector field plot using a grid described by xlim, ylim, zlim with nx, ny, and nz grid points in each direction.

Note: the vectors are represented with line, not arrow due to no implementation of :quiver3d.

F(x,y,z) = [-y, x,z]
vectorfieldplot3d(F, xlims=(-4,4), ylims=(-4,4), zlims=(0,3))

vectorfieldplot3d(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=5], [ny=5])

Create a 3 dimensional vector field plot using a grid described by xlim, ylim, zlim with nx, ny, and nz grid points in each direction.

Note: the vectors are represented with line, not arrow due to no implementation of :quiver3d.

F(x,y,z) = [-y, x,z]
vectorfieldplot3d(F, xlims=(-4,4), ylims=(-4,4), zlims=(0,3))

Plots._show(io, ::MIME"text/html", plt::Plots.Plot{Plots.PlotlyBackend})

Write a Plotly plot as its HTML body rather than as a standalone document, so the figure embeds in a rendered page — a Quarto chapter, a Jupyter cell — and stays interactive. Inherited from upstream CalculusWithJulia, which kept it behind a Plots package extension; here Plots is a hard dependency, so it lives in the package proper.

find_zero(ex, x0; kwargs...)
find_zeros(ex, a, b; kwargs...)
ZeroProblem(ex, x0)

Solve for the zeros of a symbolic expression or ~ equation, wherever Roots expects a function.

julia> using CalculusWithJuliaSquared

julia> @variables x;

julia> find_zero(x^3 - x + 1, (-2, -1))
-1.324717957244746

julia> find_zero(cos(x) ~ x, (0, 2))
0.7390851332151607

julia> find_zeros(x^2 - 1, -3, 3)
2-element Vector{Float64}:
 -1.0
  1.0

The expression must contain exactly one free variable, since nothing in the call names the one being solved for; substitute values for the others first, e.g. substitute(ex, Dict(a => 1)). An equation lhs ~ rhs is solved as lhs - rhs == 0.

Roots ships precisely this for SymPy (RootsSymPyExt) but has no Symbolics equivalent. These methods are that extension's mirror image, with Symbolics.build_function in place of lambdify.