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

Merge pull request #27 from jverzani/v0.5

Browse Source 
V0.5
This commit is contained in:
john verzani
2022-08-26 14:23:47 -04:00
committed by GitHub
parent 1b2b6cc205 4052532f38
commit 7eeea10520
21 changed files with 1751 additions and 527 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
quarto/publish.yml → .github/workflows/publish.yml vendored
View File

@@ -1,5 +1,5 @@
# This is *if* we have CI publish; currently we use `quarto render gh-pages`
on:
workflow_dispatch:
push:
branches: main
@@ -17,10 +17,18 @@ jobs:
- name: Set up Quarto
uses: quarto-dev/quarto-actions/setup@v2
- name: Install Python and Dependencies
uses: actions/setup-python@v4
with:
python-version: '3.10'
cache: 'pip'
- run: pip install jupyter
- name: Render and Publish
uses: quarto-dev/quarto-actions/publish@v2
with:
target: gh-pages
render: false
path: quarto/
env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

625
CwJ/alternatives/SciML.jmd Normal file
View File

@@ -0,0 +1,625 @@
# The SciML suite of packages
The `Julia` ecosystem advances rapidly. For much of it, the driving force is the [SciML](https://github.com/SciML) organization (Scientific Machine Learning).
In this section we describe some packages provided by this organization that could be used as alternatives to the ones utilized in these notes. Members of this organization created many packages for solving different types of differential equations, and have branched out from there. Many newer efforts of this organization have been to write uniform interfaces to other packages in the ecosystem, some of which are discussed below. We don't discuss the promise of SCIML: "Performance is considered a priority, and performance issues are considered bugs," as we don't pursue features like in-place modification, sparsity, etc. Interested readers should consult the relevant packages documentation.
The basic structure to use these packages is the "problem-algorithm-solve" interface described in [The problem-algorithm-solve interface](../ODEs/solve.html). We also discussed this interface a bit in [ODEs](../ODEs/differential_equations.html).
!!! note
These packages are in a process of rapid development and change to them is expected. These notes were written using the following versions:
```julia
pkgs = ["Symbolics", "NonlinearSolve", "Optimization", "Integrals"]
import Pkg; Pkg.status(pkgs)
```
## Symbolic math (`Symbolics`)
The `Symbolics`, `SymbolicUtils`, and `ModelingToolkit` packages are provided by this organization. These can be viewed as an alternative to `SymPy`, which is used throughout this set of notes. See the section on [Symbolics](./symbolics.html) for some additional details, the package [documentation](https://symbolics.juliasymbolics.org/stable/), or the documentation for [SymbolicsUtils](https://github.com/JuliaSymbolics/SymbolicUtils.jl).
## Solving equations
Solving one or more equations (simultaneously) is different in the linear case (where solutions are readily found -- though performance can distinguish approaches -- and the nonlinear case -- where for most situations, numeric approaches are required.
### `LinearSolve`
The `LinearSolve` package aims to generalize the solving of linear equations. For many cases these are simply represented as matrix equations of the form `Ax=b`, from which `Julia` (borrowing from MATLAB) offers the interface `A \ b` to yield `x`. There are scenarios that don't naturally fit this structure and perhaps problems where different tolerances need to be specified, and the `LinearSolve` package aims to provide a common interface to handle these scenarios. As this set of notes doesn't bump into such, this package is not described here. In the symbolic case, the `Symbolics.solve_for` function was described in [Symbolics](./symbolics.html).
### `NonlinearSolve`
The `NonlinearSolve` package can be seen as an alternative to the use of the `Roots` package in this set of notes. The package presents itself as "Fast implementations of root finding algorithms in Julia that satisfy the SciML common interface."
The package is loaded through the following command:
```julia
using NonlinearSolve
```
Unlike `Roots`, the package handles problems beyond the univariate case, as such the simplest problems have a little extra setup required.
For example, suppose we want to use this package to solve for zeros of ``f(x) = x^5 - x - 1``. We could do so a few different ways.
First, we need to define a `Julia` function representing `f`. We do so with:
```julia
f(u, p) = @. (u^5 - u - 1)
```
The function definition expects a container for the "`x`" variables and allows the passing of a container to hold parameters. We could have used the dotted operations for the power and each subtraction to allow vectorization of these basic math operations, as `u` is a container of values. The `@.` macro makes adding the "dots" quite easy, as illustrated above. It converts "every function call or operator in expr into a `dot call`."
A problem is set up with this function and an initial guess. The `@SVector` specification for the guess is for performance purposes and is provided by the `StaticArrays` package.
```julia
using StaticArrays
u0 = @SVector[1.0]
prob = NonlinearProblem(f, u0)
```
The problem is solved by calling `solve` with an appropriate method specified. Here we use Newton's method. The derivative of `f` is computed automatically.
```julia
soln = solve(prob, NewtonRaphson())
```
The basic interface for retrieving the solution from the solution object is to use indexing:
```julia
soln[]
```
----
!!! note
This interface is more performant than `Roots`, though it isn't an apples to oranges comparison as different stopping criteria are used by the two. In order to be so, we need to help out the call to `NonlinearProblem` to indicate the problem is non-mutating by adding a "`false`", as follows:
```julia
using BenchmarkTools
@btime solve(NonlinearProblem{false}(f, @SVector[1.0]), NewtonRaphson())
```
As compared to:
```julia
import Roots
import ForwardDiff
g(x) = x^5 - x - 1
gp(x) = ForwardDiff.derivative(g, x)
@btime solve(Roots.ZeroProblem((g, gp), 1.0), Roots.Newton())
```
----
This problem can also be solved using a bracketing method. The package has both `Bisection` and `Falsi` as possible methods. To use a bracketing method, the initial bracket must be specified.
```julia
u0 = (1.0, 2.0)
prob = NonlinearProblem(f, u0)
```
And
```julia
solve(prob, Bisection())
```
----
Incorporating parameters is readily done. For example to solve ``f(x) = \cos(x) - x/p`` for different values of ``p`` we might have:
```julia
f(x, p) = @. cos(x) - x/p
u0 = (0, pi/2)
p = 2
prob = NonlinearProblem(f, u0, p)
solve(prob, Bisection())
```
!!! note
The *insignificant* difference in stopping criteria used by `NonlinearSolve` and `Roots` is illustrated in this example, where the value returned by `NonlinearSolve` differs by one floating point value:
```julia
an = solve(NonlinearProblem{false}(f, u0, p), Bisection())
ar = solve(Roots.ZeroProblem(f, u0), Roots.Bisection(); p=p)
nextfloat(an[]) == ar, f(an[], p), f(ar, p)
```
----
We can solve for several parameters at once, by using an equal number of initial positions as follows:
```julia
ps = [1, 2, 3, 4]
u0 = @SVector[1, 1, 1, 1]
prob = NonlinearProblem(f, u0, ps)
solve(prob, NewtonRaphson())
```
### Higher dimensions
We solve now for a point on the surface of the following `peaks` function where the gradient is ``0``. (The gradient here will be a vector-valued function from ``R^2 `` to ``R^2.``) First we define the function:
```julia
function _peaks(x, y)
p = 3 * (1 - x)^2 * exp(-x^2 - (y + 1)^2)
p -= 10 * (x / 5 - x^3 - y^5) * exp(-x^2 - y^2)
p -= 1/3 * exp(-(x + 1)^2 - y^2)
p
end
peaks(u) = _peaks(u[1], u[2]) # pass container, take first two components
```
The gradient can be computed different ways within `Julia`, but here we use the fact that the `ForwardDiff` package is loaded by `NonlinearSolve`. Once the function is defined, the pattern is similar to above. We provide a starting point, create a problem, then solve:
```julia
∇peaks(x, p=nothing) = NonlinearSolve.ForwardDiff.gradient(peaks, x)
u0 = @SVector[1.0, 1.0]
prob = NonlinearProblem(∇peaks, u0)
u = solve(prob, NewtonRaphson())
```
We can see that this identified value is a "zero" through:
```julia; error=true
∇peaks(u.u)
```
### Using Modeling toolkit to model the non-linear problem
Nonlinear problems can also be approached symbolically using the `ModelingToolkit` package. There is one additional step necessary.
As an example, we look to solve numerically for the zeros of ``x^5-x-\alpha`` for a parameter ``\alpha``. We can describe this equation as follows:
```julia
using ModelingToolkit
@variables x
@parameters α
eq = x^5 - x - α ~ 0
```
The extra step is to specify a "`NonlinearSystem`." It is a system, as in practice one or more equations can be considered. The `NonlinearSystem`constructor handles the details where the equation, the variable, and the parameter are specified. Below this is done using vectors with just one element:
```julia
ns = NonlinearSystem([eq], [x], [α], name=:ns)
```
The `name` argument is special. The name of the object (`ns`) is assigned through `=`, but the system must also know this same name. However, the name on the left is not known when the name on the right is needed, so it is up to the user to keep them synchronized. The `@named` macro handles this behind the scenes by simply rewriting the syntax of the assignment:
```julia
@named ns = NonlinearSystem([eq], [x], [α])
```
With the system defined, we can pass this to `NonlinearProblem`, as was done with a function. The parameter is specified here, and in this case is `α => 1.0`. The initial guess is `[1.0]`:
```julia
prob = NonlinearProblem(ns, [1.0], [α => 1.0])
```
The problem is solved as before:
```julia
solve(prob, NewtonRaphson())
```
## Optimization (`Optimization.jl`)
We describe briefly the `Optimization` package which provides a common interface to *numerous* optimization packages in the `Julia` ecosystem. We discuss only the interface for `Optim.jl` defined in `OptimizationOptimJL`.
We begin with a simple example from first semester calculus:
> Among all rectangles of fixed perimeter, find the one with the *maximum* area.
If the perimeter is taken to be ``25``, the mathematical setup has a
constraint (``P=25=2x+2y``) and an objective (``A=xy``) to
maximize. In this case, the function to *maximize* is ``A(x) = x \cdot
(25-2x)/2``. This is easily done different ways, such as finding the
one critical point and identifying this as the point of maximum.
To do this last step using `Optimization` we would have.
```julia
height(x) = @. (25 - 2x)/2
A(x, p=nothing) = @.(- x * height(x))
```
The minus sign is needed here as optimization routines find *minimums*, not maximums.
To use `Optimization` we must load the package **and** the underlying backend glue code we intend to use:
```julia
using Optimization
using OptimizationOptimJL
```
Next, we define an optimization function with information on how its derivatives will be taken. The following uses `ForwardDiff`, which is a good choice in the typical calculus setting, where there are a small number of inputs (just ``1`` here.)
```julia
F = OptimizationFunction(A, Optimization.AutoForwardDiff())
x0 = [4.0]
prob = OptimizationProblem(F, x0)
```
The problem is solved through the common interface with a specified method, in this case `Newton`:
```julia
soln = solve(prob, Newton())
```
!!! note
We use `Newton` not `NewtonRaphson` as above. Both methods are similar, but they come from different uses -- for latter for solving non-linear equation(s), the former for solving optimization problems.
The solution is an object containing the identified answer and more. To get the value, use index notation:
```julia
soln[]
```
The corresponding ``y`` value and area are found by:
```julia
xstar = soln[]
height(xstar), A(xstar)
```
The `minimum` property also holds the identified minimum:
```julia
soln.minimum # compare with A(soln[], nothing)
```
The package is a wrapper around other packages. The output of the underlying package is presented in the `original` property:
```julia
soln.original
```
----
This problem can also be approached symbolically, using `ModelingToolkit`.
For example, we set up the problem with:
```julia
using ModelingToolkit
@parameters P
@variables x
y = (P - 2x)/2
Area = - x*y
```
The above should be self explanatory. To put into a form to pass to `solve` we define a "system" by specifying our objective function, the variables, and the parameters.
```julia
@named sys = OptimizationSystem(Area, [x], [P])
```
(This step is different, as before an `OptimizationFunction` was defined; we use `@named`, as above, to ensure the system has the same name as the identifier, `sys`.)
This system is passed to `OptimizationProblem` along with a specification of the initial condition (``x=4``) and the perimeter (``P=25``). A vector of pairs is used below:
```julia
prob = OptimizationProblem(sys, [x => 4.0], [P => 25.0]; grad=true, hess=true)
```
The keywords `grad=true` and `hess=true` instruct for automatic derivatives to be taken as needed. These are needed in the choice of method, `Newton`, below.
Solving this problem then follows the same pattern as before, again with `Newton` we have:
```julia
solve(prob, Newton())
```
(A derivative-free method like `NelderMead()` could be used and then the `grad` and `hess` keywords above would be unnecessary, though not harmful.)
----
The related calculus problem:
> Among all rectangles with a fixed area, find the one with *minimum* perimeter
could be similarly approached:
```julia
@parameters Area
@variables x
y = Area/x # from A = xy
P = 2x + 2y
@named sys = OptimizationSystem(P, [x], [Area])
u0 = [x => 4.0]
p = [Area => 25.0]
prob = OptimizationProblem(sys, u0, p; grad=true, hess=true)
soln = solve(prob, LBFGS())
```
We used an initial guess of ``x=4`` above. The `LBFGS` method is
a computationally efficient modification of the
Broyden-Fletcher-Goldfarb-Shanno algorithm ... It is a quasi-Newton
method that updates an approximation to the Hessian using past
approximations as well as the gradient." On this problem it performs similarly to `Newton`, though in general may be preferable for higher-dimensional problems.
### Two dimensional
Scalar functions of two input variables can have their minimum value identified in the same manner using `Optimization.jl`.
For example, consider the function
```math
f(x,y) = (x + 2y - 7)^2 + (2x + y - 5)^2
```
We wish to minimize this function.
We begin by defining a function in `Julia`:
```julia
function f(u, p)
x, y = u
(x + 2y - 7)^2 + (2x + y - 5)^2
end
```
We turn this into an optimization function by specifying how derivatives will be taken, as we will the `LBFGS` algorithm below:
```julia
ff = OptimizationFunction(f, Optimization.AutoForwardDiff())
```
We will begin our search at the origin. We have to be mindful to use floating point numbers here:
```julia
u0 = [0.0, 0.0] # or zeros(2)
```
```julia
prob = OptimizationProblem(ff, u0)
```
Finally, we solve the values:
```julia
solve(prob, LBFGS())
```
The value of ``(1, 3)`` agrees with the contour graph, as it is a point in the interior of the contour for the smallest values displayed.
```julia
using Plots
xs = range(0, 2, length=100)
ys = range(2, 4, length=100)
contour(xs, ys, (x,y) -> f((x,y), nothing))
```
We could also use a *derivative-free* method, and skip a step:
```julia
prob = OptimizationProblem(f, u0) # skip making an OptimizationFunction
solve(prob, NelderMead())
```
## Integration (`Integrals.jl`)
The `Integrals` package provides a common interface to different numeric integration packages in the `Julia` ecosystem. For example, `QuadGK` and `HCubature`. The value of this interface, over those two packages, is its non-differentiated access to other packages, which for some uses may be more performant.
The package follows the same `problem-algorithm-solve` interface, as already seen.
The interface is designed for ``1``-and-higher dimensional integrals.
The package is loaded with
```julia
using Integrals
```
For a simple definite integral, such as ``\int_0^\pi \sin(x)dx``, we have:
```julia
f(x, p) = sin(x)
prob = IntegralProblem(f, 0.0, pi)
soln = solve(prob, QuadGKJL())
```
To get access to the answer, we can use indexing notation:
```julia
soln[]
```
Comparing to just using `QuadGK`, the same definite integral would be
estimated with:
```julia
using QuadGK
quadgk(sin, 0, pi)
```
The estimated upper bound on the error from `QuadGK`, is available through the `resid` property on the `Integrals` output:
```julia
soln.resid
```
The `Integrals` solution is a bit more verbose, but it is more flexible. For example, the `HCubature` package provides a similar means to compute ``n``- dimensional integrals. For this problem, the modifications would be:
```julia
f(x, p) = sin.(x)
prob = IntegralProblem(f, [0.0], [pi])
soln = solve(prob, HCubatureJL())
```
```julia
soln[]
```
The estimated maximum error is also given by `resid`:
```julia
soln.resid
```
----
As well, suppose we wanted to parameterize our function and then differentiate.
Consider ``d/dp \int_0^\pi \sin(px) dx``. We can do this integral directly to get
```math
\begin{align*}
\frac{d}{dp} \int_0^\pi \sin(px)dx
&= \frac{d}{dp}\left( \frac{-1}{p} \cos(px)\Big\rvert_0^\pi\right)\\
&= \frac{d}{dp}\left( -\frac{\cos(p\cdot\pi)-1}{p}\right)\\
&= \frac{\cos(p\cdot \pi) - 1)}{p^2} + \frac{\pi\cdot\sin(p\cdot\pi)}{p}
\end{align*}
```
Using `Integrals` with `QuadGK` we have:
```julia
f(x, p) = sin(p*x)
function ∫sinpx(p)
prob = IntegralProblem(f, 0.0, pi, p)
solve(prob, QuadGKJL())
end
```
We can compute values at both ``p=1`` and ``p=2``:
```julia
∫sinpx(1), ∫sinpx(2)
```
To find the derivative in ``p`` , we have:
```julia
ForwardDiff.derivative(∫sinpx, 1), ForwardDiff.derivative(∫sinpx, 2)
```
(In `QuadGK`, the following can be differentiated `∫sinpx(p) = quadgk(x -> sin(p*x), 0, pi)[1]` as well. `Integrals` gives a consistent interface.
### Higher dimension integrals
The power of a common interface is the ability to swap backends and the uniformity for different dimensions. Here we discuss integrals of scalar-valued and vector-valued functions.
#### ``f: R^n \rightarrow R``
The area under a surface generated by ``z=f(x,y)`` over a rectangular region ``[a,b]\times[c,d]`` can be readily computed. The two coding implementations require ``f`` to be expressed as a function of a vector--*and* a parameter--and the interval to be expressed using two vectors, one for the left endpoints (`[a,c]`) and on for the right endpoints (`[b,d]`).
For example, the area under the function ``f(x,y) = 1 + x^2 + 2y^2`` over ``[-1/2, 1/2] \times [-1,1]`` is computed by:
```julia
f(x, y) = 1 + x^2 + 2y^2 # match math
fxp(x, p) = f(x[1], x[2]) # prepare for IntegralProblem
ls = [-1/2, -1] # left endpoints
rs = [1/2, 1] # right endpoints
prob = IntegralProblem(fxp, ls, rs)
soln = solve(prob, HCubatureJL())
```
Of course, we could have directly defined the function (`fxp`) using indexing of the `x` variable.
----
For non-rectangular domains a change of variable is required.
For example, an integral to assist in finding the volume of a sphere might be
```math
V = 2 \iint_R \sqrt{\rho^2 - x^2 - y^2} dx dy
```
where ``R`` is the disc of radius ``\rho`` in the ``x-y`` plane.
The usual approach is to change to polar-coordinates and write this integral as
```math
V = \int_0^{2\pi}\int_0^\rho \sqrt{\rho^2 - r^2} r dr d\theta
```
the latter being an integral over a rectangular domain.
To compute this transformed integral, we might have:
```julia
function vol_sphere(ρ)
f(rθ, p) = sqrt(ρ^2 - rθ[1]^2) * rθ[1]
ls = [0,0]
rs = [ρ, 2pi]
prob = IntegralProblem(f, ls, rs)
solve(prob, HCubatureJL())
end
vol_sphere(2)
```
If it is possible to express the region to integrate as ``G(R)`` where ``R`` is a rectangular region, then the change of variables formula,
```math
\iint_{G(R)} f(x) dA = \iint_R (f\circ G)(u) |det(J_G(u)| dU
```
turns the integral into the non-rectangular domain ``G(R)`` into one over the rectangular domain ``R``. The key is to *identify* ``G`` and to compute the Jacobian. The latter is simply accomplished with `ForwardDiff.jacobian`.
For an example, we find the moment of inertia about the axis of the unit square tilted counter-clockwise an angle ``0 \leq \alpha \leq \pi/2``.
The counter clockwise rotation of a unit square by angle ``\alpha`` is described by
```math
G(u, v) = \langle \cos(\alpha)\cdot u - \sin(\alpha)\cdot v, \sin(\alpha)\cdot u, +\cos(\alpha)\cdot v \rangle
```
So we have ``\iint_{G(R)} x^2 dA`` is computed by the following with ``\alpha=\pi/4``:
```julia
import LinearAlgebra: det
𝑓(uv) = uv[1]^2
function G(uv)
α = pi/4 # could be made a parameter
u,v = uv
[cos(α)*u - sin(α)*v, sin(α)*u + cos(α)*v]
end
f(u, p) = (𝑓∘G)(u) * det(ForwardDiff.jacobian(G, u))
prob = IntegralProblem(f, [0,0], [1,1])
solve(prob, HCubatureJL())
```
#### ``f: R^n \rightarrow R^m``
The `Integrals` package provides an interface for vector-valued functions. By default, the number of dimensions in the output is assumed to be ``1``, but the `nout` argument can adjust that.
Let ``f`` be vector valued with components ``f_1, f_2, \dots, f_m``, then the output below is the vector with components ``\iint_R f_1 dV, \iint_R f_2 dV, \dots, \iint_R f_m dV``.
For a trivial example, we have:
```julia
f(x, p) = [x[1], x[2]^2]
prob = IntegralProblem(f, [0,0],[3,4], nout=2)
solve(prob, HCubatureJL())
```

10
CwJ/alternatives/sciml_interface.__jmd__
View File

@@ -1,10 +0,0 @@
# seems like with
* quadrature (renamed)
* nonlinsolve
* GalacticOptim (renamed)
* symbolic-numeric integration
* symbolics.jl
...
This should be mentioned

230
CwJ/alternatives/symbolics.jmd
View File

@@ -4,7 +4,7 @@ There are a few options in `Julia` for symbolic math, for example, the `SymPy` p
## About
The `Symbolics` package bill itself as a "fast and modern Computer Algebra System (CAS) for a fast and modern programming language." This package relies on the `SymbolicUtils` package and is built upon by the `ModelingToolkit` package, which we don't describe here.
The `Symbolics` package bills itself as a "fast and modern Computer Algebra System (CAS) for a fast and modern programming language." This package relies on the `SymbolicUtils` package and is built upon by the `ModelingToolkit` package, which is only briefly touched on here.
We begin by loading the `Symbolics` package which when loaded re-exports the `SymbolicUtils` package.
@@ -15,7 +15,7 @@ using Symbolics
## Symbolic variables
Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two package discussed here.
Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two packages discussed here.
`SymbolicUtils` creates variables which carry `Julia` type information (e.g. `Int`, `Float64`, ...). This type information carries through operations involving these variables. Symbolic variables can be created with the `@syms` macro. For example
@@ -23,9 +23,9 @@ Symbolic math at its core involves symbolic variables, which essentially defer e
@syms x y::Int f(x::Real)::Real
```
This creates `x` a symbolic value with type `Number`, `y` a symbolic variable holding integer values, and `f` a symbolic function of a single real variable outputing a real variable.
This creates `x` a symbolic value with symbolic type `Number`, `y` a symbolic variable holding integer values, and `f` a symbolic function of a single real variable outputting a real variable.
The `symtype` function reveals the underlying type:
The non-exported `symtype` function reveals the underlying type:
```julia
import Symbolics.SymbolicUtils: symtype
@@ -33,14 +33,14 @@ import Symbolics.SymbolicUtils: symtype
symtype(x), symtype(y)
```
For `y`, the symbolic type being real does not imply the `y` has a subtype of `Real`:
For `y`, the symbolic type being real does not imply the type of `y` is a subtype of `Real`:
```julia
isa(y, Real)
```
We see that the function `f`, when called with `y` would return a value of (symbolic) type `Real`:
We see that the function `f` when called with `y` would return a value of (symbolic) type `Real`:
```julia
f(y) |> symtype
@@ -55,8 +55,9 @@ f(x)
!!! note
The `SymPy` package also has an `@syms` macro to create variables. Though they names agree, they do different things. Using both packages together would require qualifying many shared method names. For `SymbolicUtils`, the `@syms` macro uses `Julia` types to parameterize the variables. In `SymPy` it is possible to specify *assumptions* on the variables, but that is different and not useful for dispatch without some extra effort.
The `SymPy` package also has an `@syms` macro to create variables. Though their names agree, they do different things. Using both packages together would require qualifying many shared method names. For `SymbolicUtils`, the `@syms` macro uses `Julia` types to parameterize the variables. In `SymPy` it is possible to specify *assumptions* on the variables, but that is different and not useful for dispatch without some extra effort.
### Variables in Symbolics
For `Symbolics`, symbolic variables are created using a wrapper around an underlying `SymbolicUtils` object. This wrapper, `Num`, is a subtype of `Real` (the underlying `SymbolicUtils` object may have symbolic type `Real`, but it won't be a subtype of `Real`.)
@@ -79,6 +80,41 @@ typeof(x), symtype(x), typeof(Symbolics.value(x))
```
(The `value` method unwraps the `Num` wrapper.)
### Variables in ModelingToolkit
The `ModelingToolkit` package has a slightly different declaration for variables, and is described next. First the package is loaded
```julia
using ModelingToolkit
```
`ModelingToolkit` re-exports all of the `Symbolics` package when loaded.
The role of `ModelingToolkit` is that "it allows for users to give a high-level description of a model for symbolic preprocessing to analyze and enhance the model." This symbolic description allows for variables to be identified as "parameters" or "variables". For example, to parameterize a quadratic equation:
```julia
@parameters a b c
@variables x
y = a*x^2 + b*x + c
```
The numeric solution of the quadratic equation (solving for ``y=0``) would involved specifying values for the parameters and then numerically solving. This separation of parameters and variables is similar to the `f(x, p)` pattern of function definition.
The typical usage is multi-variable. This example is from the package's documentationto describe a differential equation:
```julia
@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)
```
The `D` will be described a bit later, but it formally specifies a derivative in the `t` variable. The `x(t)`, `y(t)`, `z(t)` are symbolic functions of `t`, so expressions like `D(x)` below mean the time derivative of an unknown function `x`. Here are how the Lorenz equation equations are specified:
```julia
eqs = [D(x) ~ σ * (y - x),
D(y) ~ x * (ρ - z) - y,
D(z) ~ x * y - β * z]
```
## Symbolic expressions
@@ -110,7 +146,7 @@ typeof(sin(x)), typeof(Symbolics.value(sin(x)))
### Tree structure to expressions
The `TermInterface` package is used by `SymbolicUtils` to explore the tree structdure of an expression. The main methods are (cf. [symbolicutils.jl](https://symbolicutils.juliasymbolics.org/#expression_interface)):
The `TermInterface` package is used by `SymbolicUtils` to explore the tree structure of an expression. The main methods are (cf. [SymbolicUtils.jl](https://symbolicutils.juliasymbolics.org/#expression_interface)):
* `istree(ex)`: `true` if `ex` is not a *leaf* node (like a symbol or numeric literal)
* `operation(ex)`: the function being called (if `istree` returns `true`)
@@ -184,7 +220,7 @@ free_symbols(ex)
### Substitute
The `substitute` command is used to replace values with other values. For examples:
The `substitute` command is used to replace values with other values. For example:
```julia
@variables x y z
@@ -217,20 +253,20 @@ substitute(ex, x=>π), substitute(ex, x=>π, fold=false)
Algebraic operations with symbolic values can involve an exponentially increasing number of terms. As such, some simplification rules are applied after an operation to reduce the complexity of the computed value.
For example, `0+x` should simplify to `x`, `0+x`, `x^0`, or `x^1` should simplify, to some natural answer.
For example, `0+x` should simplify to `x`, as well `1*x`, `x^0`, or `x^1` should each simplify, to some natural answer.
`SymbolicUtils` also [simplifies](https://symbolicutils.juliasymbolics.org/#simplification) several other expressions, including:
* `-x` becomes `(-1)*x`
* `x * x` becomes `x^2` (and `x^n` if more terms). Meaning this expression is represented as a power, not a product
* `x + x` becomes `2*x` (and `n*x` if more terms). Similarly, this represented as a product, not a sum.
* `p/q * x` becomes `(p*x)/q)`, similarly `p/q * x/y` becomes `(p*x)/(q*y)`
* `p/q * x` becomes `(p*x)/q)`, similarly `p/q * x/y` becomes `(p*x)/(q*y)`. (Division wraps multiplication.)
In `SymbolicUtils`, this *rewriting* is accomplished by means of *rewrite rules*. The package makes it easy to apply user-written rewrite rules.
### Rewriting
Many algebraic simplifications are done by the `simplify` command. For example, the basic trignometric identities are applied:
Many algebraic simplifications are done by the `simplify` command. For example, the basic trigonometric identities are applied:
```julia
@variables x
@@ -255,10 +291,12 @@ ex1 = substitute(ex, x => sin(x + y + z))
ex1 |> Symbolics.value |> r |> Num
```
Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trignometry and exponential simplifications:
Rewrite rules when applied return the rewritten expression, if there is a match, or `nothing`.
Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trigonometry and exponential simplifications:
```julia
r = @rule(exp(~x)^(~y) => exp(~x * ~y))
r = @rule(exp(~x)^(~y) => exp(~x * ~y)) # (e^x)^y -> e^(x*y)
ex = exp(-x+z)^y
ex, ex |> Symbolics.value |> r |> Num
```
@@ -289,7 +327,7 @@ ex = exp(-x + z)^y
Symbolics.toexpr(ex)
```
This output shows an internal representation of the steps for computing the value `ex` given different inputs.
This output shows an internal representation of the steps for computing the value `ex` given different inputs. (The number `(-1)` multiplies `x`, this is added to `z` and the result passed to `exp`. That values is then used as the base for `^` with exponent `y`.)
Such `Julia` expressions are one step away from building `Julia` functions for evaluating symbolic expressions fast (though with some technical details about "world age" to be reckoned with). The `build_function` function with the argument `expression=Val(false)` will compile a `Julia` function:
@@ -311,7 +349,7 @@ However, `build_function` will be **significantly** more performant, which when
The above, through passing ``3`` variables after the expression, creates a function of ``3`` variables. Functions of a vector of inputs can also be created, just by expressing the variables in that manner:
```juila
```julia
h1 = build_function(ex, [x, y, z]; expression=Val(false))
h1([1, 2, 3]) # not h1(1,2,3)
```
@@ -330,7 +368,7 @@ Roots.find_zero(λ, (1, 2))
### Plotting
Using `Plots`, the plotting symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via `build_function`.
Using `Plots`, the plotting of symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via `build_function`.
For example,
@@ -376,7 +414,7 @@ For example
d, r = polynomial_coeffs(a*x^2 + b*x + c, (x,))
```
The first term output is dictionary who's keys are the monomials and who's values are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just the constant `c`.
The first term output is dictionary with keys which are the monomials and with values which are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just the constant `c`.
The expression can then be reconstructed through
@@ -384,7 +422,7 @@ The expression can then be reconstructed through
r + sum(v*k for (k,v) ∈ d)
```
The above has `a,b,c` as parameters and `x` as the symbol. This separation is designated by passing the desired polynomials symbols to `polynomial_coeff` as an iterable. (Above as a ``1``-element tuple.)
The above has `a,b,c` as parameters and `x` as the symbol. This separation is designated by passing the desired polynomial symbols to `polynomial_coeff` as an iterable. (Above as a ``1``-element tuple.)
More complicated polynomials can be similarly decomposed:
@@ -464,6 +502,62 @@ m,n = degree.(nd(ex))
m > n ? "limit is infinite" : m < n ? "limit is 0" : "limit is a constant"
```
### Vectors and matrices
Symbolic vectors and matrices can be created with a specified size:
```julia
@variables v[1:3] M[1:2, 1:3] N[1:3, 1:3]
```
Computations, like finding the determinant below, are lazy unless the values are `collect`ed:
```julia
using LinearAlgebra
det(N)
```
```julia
det(collect(N))
```
Similarly, with `norm`:
```julia
norm(v)
```
and
```julia
norm(collect(v))
```
Matrix multiplication is also deferred, but the size compatability of the matrices and vectors is considered early:
```julia
M*N, N*N, M*v
```
This errors, as the matrix dimensions are not compatible for multiplication:
```julia; error=true
N*M
```
Similarly, linear solutions can be symbolically specified:
```julia
@variables R[1:2, 1:2] b[1:2]
R \ b
```
```julia
collect(R \ b)
```
### Algebraically solving equations
The `~` operator creates a symbolic equation. For example
@@ -476,15 +570,33 @@ x^5 - x ~ 1
or
```julia
ex = [5x + 2y, 6x + 3y] .~ [1, 2]
eqs = [5x + 2y, 6x + 3y] .~ [1, 2]
```
The `Symbolics.solve_for` function can solve *linear* equations. For example,
```julia
Symbolics.solve_for(ex, [x, y])
Symbolics.solve_for(eqs, [x, y])
```
The coefficients can be symbolic. Two examples could be:
```julia
@variables m b x y
eq = y ~ m*x + b
Symbolics.solve_for(eq, x)
```
```julia
@variables a11 a12 a22 x y b1 b2
R,X,b = [a11 a12; 0 a22], [x; y], [b1, b2]
eqs = R*X .~ b
```
```julia
Symbolics.solve_for(eqs, [x,y])
```
### Limits
@@ -496,24 +608,41 @@ As of writing, there is no extra functionality provided by `Symbolics` for compu
```julia
@variables a b c x
ex = a*x^2 + b*x + c
Symbolics.derivative(ex, x)
y = a*x^2 + b*x + c
yp = Symbolics.derivative(y, x)
```
The computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:
Or to find a critical point:
```julia
Symbolics.solve_for(yp ~ 0, x) # linear equation to solve
```
The derivative computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:
```julia
D = Differential(x)
D(ex)
D(y)
```
and then
```julia
expand_derivatives(D(ex))
expand_derivatives(D(y))
```
The differentials can be multiplied to create operators for taking higher-order derivatives:
Using `Differential`, differential equations can be specified. An example was given in [ODEs](../ODEs/differential_equations.html), using `ModelingToolkit`.
Higher order derivatives can be done through composition:
```julia
D(D(y)) |> expand_derivatives
```
Differentials can also be multiplied to create operators for taking higher-order derivatives:
```julia
@variables x y
@@ -529,7 +658,7 @@ In addition to `Symbolics.derivative` there are also the helper functions, such
Symbolics.hessian(ex, [x,y])
```
The `gradient` function is also available
The `gradient` function is also defined
```julia
@variables x y z
@@ -537,12 +666,12 @@ ex = x^2 - 2x*y + z*y
Symbolics.gradient(ex, [x, y, z])
```
The `jacobian` takes an array of expressions:
The `jacobian` function takes an array of expressions:
```julia
@variables x y
exs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(exs, [x,y])
eqs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(eqs, [x,y])
```
@@ -553,12 +682,14 @@ The `SymbolicNumericIntegration` package provides a means to integrate *univaria
Symbolic integration can be approached in different ways. SymPy implements part of the Risch algorithm in addition to other algorithms. Rules-based algorithms could also be implemented.
For example, here is a simple rule that could be used to integrate a single integral
For a trivial example, here is a rule that could be used to integrate a single integral
```julia
is_var(x) = (xs = Symbolics.get_variables(x); length(xs) == 1 && xs[1] === x)
@syms x ∫(x)
is_var(x) = (xs = Symbolics.get_variables(x); length(xs) == 1 && xs[1] === x)
r = @rule ∫(~x::is_var) => x^2/2
r(∫(x))
```
@@ -567,19 +698,34 @@ The `SymbolicNumericIntegration` package includes many more predicates for doing
If ``f(x)`` is to be integrated, a set of *candidate* answers is generated. The following is **proposed** as an answer: ``\sum q_i \Theta_i(x)``. Differentiating the proposed answer leads to a *linear system of equations* that can be solved.
The example in the [paper](https://arxiv.org/pdf/2201.12468v2.pdf) describing the method is with ``f(x) = x \sin(x)`` and the candidate thetas are ``{x, \sin(x), \cos(x), x\sin(x), x\cos(x)}`` so that we propose:
The example in the [paper](https://arxiv.org/pdf/2201.12468v2.pdf) describing the method is with ``f(x) = x \sin(x)`` and the candidate thetas are ``{x, \sin(x), \cos(x), x\sin(x), x\cos(x)}`` so that the propose answer is:
```math
\int f(x) dx = q_1 x + q_2 \sin(x) + q_3 \cos(x) + q_4 x \sin(x) + q_4 x \cos(x)
```
Differentiating both sides, yields a term ``x\sin(x)`` on the left, and equating coefficients gives:
We differentiate the right hand side:
```math
q_1 = q_4 = 0,\quad q_5 = -1, \quad q_4 - q_3 = q_2 - q_5 = 0
```julia
@variables q[1:5] x
ΣqᵢΘᵢ = dot(collect(q), (x, sin(x), cos(x), x*sin(x), x*cos(x)))
simplify(Symbolics.derivative(ΣqᵢΘᵢ, x))
```
which can be solved with ``q_5=-1``, ``q_2=1``, and the other coefficients being ``0``. That is ``\int f(x) dx = 1 \sin(x) + (-1) x\cos(x)``.
This must match ``x\sin(x)`` so we have by
equating coefficients of the respective terms:
```math
q_2 + q_5 = 0, \quad q_4 = 0, \quad q_1 = 0, \quad q_3 = 0, \quad q_5 = -1
```
That is ``q_2=1``, ``q_5=-1``, and the other coefficients are ``0``, giving
an answer computed with:
```julia
d = Dict(q[i] => v for (i,v) ∈ enumerate((0,1,0,0,-1)))
substitute(ΣqᵢΘᵢ, d)
```
The package provides an algorithm for the creation of candidates and the means to solve when possible. The `integrate` function is the main entry point. It returns three values: `solved`, `unsolved`, and `err`. The `unsolved` is the part of the integrand which can not be solved through this package. It is `0` for a given problem when `integrate` is successful in identifying an antiderivative, in which case `solved` is the answer. The value of `err` is a bound on the numerical error introduced by the algorithm.
@@ -592,7 +738,7 @@ using SymbolicNumericIntegration
integrate(x * sin(x))
```
The second term is `0`, as this has an identified antiderivative.
The second term is `0`, as this integrand has an identified antiderivative.
```julia
integrate(exp(x^2) + sin(x))
@@ -612,6 +758,8 @@ The derivative of `u` matches up to some numeric tolerance:
Symbolics.derivative(u, x) - sin(x)^5
```
----
The integration of rational functions (ratios of polynomials) can be done algorithmically, provided the underlying factorizations can be identified. The `SymbolicNumericIntegration` package has a function `factor_rational` that can identify factorizations.
```julia
@@ -641,13 +789,13 @@ u = 1 / expand((x^2+1)*(x-2)^2)
v = factor_rational(u)
```
As such, the integrals have numeric differences:
As such, the integrals have numeric differences from their mathematical counterparts:
```julia
a,b,c = integrate(u)
```
We can see a bit of why through the following which needs a tolerance set to identify the rational numbers correctly:
We can see a bit of how much through the following, which needs a tolerance set to identify the rational numbers of the mathematical factorization correctly:
```julia
cs = [first(arguments(term)) for term ∈ arguments(a)] # pick off coefficients

4
CwJ/precalc/transformations.jmd
View File

@@ -359,8 +359,8 @@ earth around the sun, is not accurate enough for precise work, but it does help
##### Example: a growth model in fisheries
The von Bertanlaffy growth
[equation](http://www.fao.org/docrep/W5449e/w5449e05.htm) is
The von Bertalanffy growth
[equation](https://en.wikipedia.org/wiki/Von_Bertalanffy_function) is
$L(t) =L_\infty \cdot (1 - e^{k\cdot(t-t_0)})$. This family of functions can
be viewed as a transformation of the exponential function $f(t)=e^t$.
Part is a scaling and shifting (the $e^{k \cdot (t - t_0)}$)

2
Project.toml
View File

@@ -1,7 +1,7 @@
name = "CalculusWithJuliaNotes"
uuid = "8cd3c377-0a30-4ec5-b85a-75291d749efe"
authors = ["jverzani <jverzani@gmail.com> and contributors"]
version = "0.1.2"
version = "0.1.3"
[compat]
julia = "1"

6
quarto/_quarto.yml
View File

@@ -1,4 +1,4 @@
version: 0.4
version: 0.5
project:
type: book
@@ -108,10 +108,10 @@ book:
- integral_vector_calculus/stokes_theorem.qmd
- integral_vector_calculus/review.qmd
- part: "Alternative packages"
- part: alternatives.qmd
chapters:
- alternatives/symbolics.qmd
# - alternatives/sciML.qmd
- alternatives/SciML.qmd
# - alternatives/interval_arithmetic.qmd
- alternatives/plotly_plotting.qmd
- alternatives/makie_plotting.qmd

13
quarto/alternatives.qmd Normal file
View File

@@ -0,0 +1,13 @@
# Alternative packages
These notes use a particular selection of packages. This selection could have been different. For example:
* The symbolic math is provided by `SymPy`. [Symbolics](./alternatives/Symbolics.html) (along with `SymbolicUtils` and `ModelingToolkit`) provides an alternative.
* 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 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.
* 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.

11
quarto/alternatives/README
View File

@@ -1,11 +0,0 @@
# Alternatives
There are many ways to do related things in `Julia`. This directory holds alternatives to the some choices made within these notes:
## Symbolics
* needs writing
## Makie
* needs updating

734
quarto/alternatives/SciML.qmd Normal file
View File

@@ -0,0 +1,734 @@
# The SciML suite of packages
The `Julia` ecosystem advances rapidly. For much of it, the driving force is the [SciML](https://github.com/SciML) organization (Scientific Machine Learning).
In this section we describe some packages provided by this organization that could be used as alternatives to the ones utilized in these notes. Members of this organization created many packages for solving different types of differential equations, and have branched out from there. Many newer efforts of this organization have been to write uniform interfaces to other packages in the ecosystem, some of which are discussed below. We don't discuss the promise of SCIML: "Performance is considered a priority, and performance issues are considered bugs," as we don't pursue features like in-place modification, sparsity, etc. Interested readers should consult the relevant packages documentation.
The basic structure to use these packages is the "problem-algorithm-solve" interface described in [The problem-algorithm-solve interface](../ODEs/solve.html). We also discussed this interface a bit in [ODEs](../ODEs/differential_equations.html).
:::{.callout-note}
## Note
These packages are in a process of rapid development and change to them is expected. These notes were written using the following versions:
:::
```{julia}
pkgs = ["Symbolics", "NonlinearSolve", "Optimization", "Integrals"]
import Pkg; Pkg.status(pkgs)
```
## Symbolic math (`Symbolics`)
The `Symbolics`, `SymbolicUtils`, and `ModelingToolkit` packages are provided by this organization. These can be viewed as an alternative to `SymPy`, which is used throughout this set of notes. See the section on [Symbolics](./symbolics.html) for some additional details, the package [documentation](https://symbolics.juliasymbolics.org/stable/), or the documentation for [SymbolicsUtils](https://github.com/JuliaSymbolics/SymbolicUtils.jl).
## Solving equations
Solving one or more equations (simultaneously) is different in the linear case (where solutions are readily found – though performance can distinguish approaches – and the nonlinear case – where for most situations, numeric approaches are required.
### `LinearSolve`
The `LinearSolve` package aims to generalize the solving of linear equations. For many cases these are simply represented as matrix equations of the form `Ax=b`, from which `Julia` (borrowing from MATLAB) offers the interface `A \ b` to yield `x`. There are scenarios that don't naturally fit this structure and perhaps problems where different tolerances need to be specified, and the `LinearSolve` package aims to provide a common interface to handle these scenarios. As this set of notes doesn't bump into such, this package is not described here. In the symbolic case, the `Symbolics.solve_for` function was described in [Symbolics](./symbolics.html).
### `NonlinearSolve`
The `NonlinearSolve` package can be seen as an alternative to the use of the `Roots` package in this set of notes. The package presents itself as "Fast implementations of root finding algorithms in Julia that satisfy the SciML common interface."
The package is loaded through the following command:
```{julia}
using NonlinearSolve
```
Unlike `Roots`, the package handles problems beyond the univariate case, as such the simplest problems have a little extra setup required.
For example, suppose we want to use this package to solve for zeros of $f(x) = x^5 - x - 1$. We could do so a few different ways.
First, we need to define a `Julia` function representing `f`. We do so with:
```{julia}
f(u, p) = @. (u^5 - u - 1)
```
The function definition expects a container for the "`x`" variables and allows the passing of a container to hold parameters. We could have used the dotted operations for the power and each subtraction to allow vectorization of these basic math operations, as `u` is a container of values. The `@.` macro makes adding the "dots" quite easy, as illustrated above. It converts "every function call or operator in expr into a `dot call`."
A problem is set up with this function and an initial guess. The `@SVector` specification for the guess is for performance purposes and is provided by the `StaticArrays` package.
```{julia}
using StaticArrays
u0 = @SVector[1.0]
prob = NonlinearProblem(f, u0)
```
The problem is solved by calling `solve` with an appropriate method specified. Here we use Newton's method. The derivative of `f` is computed automatically.
```{julia}
soln = solve(prob, NewtonRaphson())
```
The basic interface for retrieving the solution from the solution object is to use indexing:
```{julia}
soln[]
```
---
:::{.callout-note}
## Note
This interface is more performant than `Roots`, though it isn't an apples to oranges comparison as different stopping criteria are used by the two. In order to be so, we need to help out the call to `NonlinearProblem` to indicate the problem is non-mutating by adding a "`false`", as follows:
:::
```{julia}
using BenchmarkTools
@btime solve(NonlinearProblem{false}(f, @SVector[1.0]), NewtonRaphson())
```
As compared to:
```{julia}
import Roots
import ForwardDiff
g(x) = x^5 - x - 1
gp(x) = ForwardDiff.derivative(g, x)
@btime solve(Roots.ZeroProblem((g, gp), 1.0), Roots.Newton())
```
---
This problem can also be solved using a bracketing method. The package has both `Bisection` and `Falsi` as possible methods. To use a bracketing method, the initial bracket must be specified.
```{julia}
u0 = (1.0, 2.0)
prob = NonlinearProblem(f, u0)
```
And
```{julia}
solve(prob, Bisection())
```
---
Incorporating parameters is readily done. For example to solve $f(x) = \cos(x) - x/p$ for different values of $p$ we might have:
```{julia}
f(x, p) = @. cos(x) - x/p
u0 = (0, pi/2)
p = 2
prob = NonlinearProblem(f, u0, p)
solve(prob, Bisection())
```
:::{.callout-note}
## Note
The *insignificant* difference in stopping criteria used by `NonlinearSolve` and `Roots` is illustrated in this example, where the value returned by `NonlinearSolve` differs by one floating point value:
:::
```{julia}
an = solve(NonlinearProblem{false}(f, u0, p), Bisection())
ar = solve(Roots.ZeroProblem(f, u0), Roots.Bisection(); p=p)
nextfloat(an[]) == ar, f(an[], p), f(ar, p)
```
---
We can solve for several parameters at once, by using an equal number of initial positions as follows:
```{julia}
ps = [1, 2, 3, 4]
u0 = @SVector[1, 1, 1, 1]
prob = NonlinearProblem(f, u0, ps)
solve(prob, NewtonRaphson())
```
### Higher dimensions
We solve now for a point on the surface of the following `peaks` function where the gradient is $0$. (The gradient here will be a vector-valued function from $R^2$ to $R^2.$) First we define the function:
```{julia}
function _peaks(x, y)
p = 3 * (1 - x)^2 * exp(-x^2 - (y + 1)^2)
p -= 10 * (x / 5 - x^3 - y^5) * exp(-x^2 - y^2)
p -= 1/3 * exp(-(x + 1)^2 - y^2)
p
end
peaks(u) = _peaks(u[1], u[2]) # pass container, take first two components
```
The gradient can be computed different ways within `Julia`, but here we use the fact that the `ForwardDiff` package is loaded by `NonlinearSolve`. Once the function is defined, the pattern is similar to above. We provide a starting point, create a problem, then solve:
```{julia}
∇peaks(x, p=nothing) = NonlinearSolve.ForwardDiff.gradient(peaks, x)
u0 = @SVector[1.0, 1.0]
prob = NonlinearProblem(∇peaks, u0)
u = solve(prob, NewtonRaphson())
```
We can see that this identified value is a "zero" through:
```{julia}
#| error: true
∇peaks(u.u)
```
### Using Modeling toolkit to model the non-linear problem
Nonlinear problems can also be approached symbolically using the `ModelingToolkit` package. There is one additional step necessary.
As an example, we look to solve numerically for the zeros of $x^5-x-\alpha$ for a parameter $\alpha$. We can describe this equation as follows:
```{julia}
using ModelingToolkit
@variables x
@parameters α
eq = x^5 - x - α ~ 0
```
The extra step is to specify a "`NonlinearSystem`." It is a system, as in practice one or more equations can be considered. The `NonlinearSystem`constructor handles the details where the equation, the variable, and the parameter are specified. Below this is done using vectors with just one element:
```{julia}
ns = NonlinearSystem([eq], [x], [α], name=:ns)
```
The `name` argument is special. The name of the object (`ns`) is assigned through `=`, but the system must also know this same name. However, the name on the left is not known when the name on the right is needed, so it is up to the user to keep them synchronized. The `@named` macro handles this behind the scenes by simply rewriting the syntax of the assignment:
```{julia}
@named ns = NonlinearSystem([eq], [x], [α])
```
With the system defined, we can pass this to `NonlinearProblem`, as was done with a function. The parameter is specified here, and in this case is `α => 1.0`. The initial guess is `[1.0]`:
```{julia}
prob = NonlinearProblem(ns, [1.0], [α => 1.0])
```
The problem is solved as before:
```{julia}
solve(prob, NewtonRaphson())
```
## Optimization (`Optimization.jl`)
We describe briefly the `Optimization` package which provides a common interface to *numerous* optimization packages in the `Julia` ecosystem. We discuss only the interface for `Optim.jl` defined in `OptimizationOptimJL`.
We begin with a simple example from first semester calculus:
> Among all rectangles of fixed perimeter, find the one with the *maximum* area.
If the perimeter is taken to be $25$, the mathematical setup has a constraint ($P=25=2x+2y$) and an objective ($A=xy$) to maximize. In this case, the function to *maximize* is $A(x) = x \cdot (25-2x)/2$. This is easily done different ways, such as finding the one critical point and identifying this as the point of maximum.
To do this last step using `Optimization` we would have.
```{julia}
height(x) = @. (25 - 2x)/2
A(x, p=nothing) = @.(- x * height(x))
```
The minus sign is needed here as optimization routines find *minimums*, not maximums.
To use `Optimization` we must load the package **and** the underlying backend glue code we intend to use:
```{julia}
using Optimization
using OptimizationOptimJL
```
Next, we define an optimization function with information on how its derivatives will be taken. The following uses `ForwardDiff`, which is a good choice in the typical calculus setting, where there are a small number of inputs (just $1$ here.)
```{julia}
F = OptimizationFunction(A, Optimization.AutoForwardDiff())
x0 = [4.0]
prob = OptimizationProblem(F, x0)
```
The problem is solved through the common interface with a specified method, in this case `Newton`:
```{julia}
soln = solve(prob, Newton())
```
:::{.callout-note}
## Note
We use `Newton` not `NewtonRaphson` as above. Both methods are similar, but they come from different uses – for latter for solving non-linear equation(s), the former for solving optimization problems.
:::
The solution is an object containing the identified answer and more. To get the value, use index notation:
```{julia}
soln[]
```
The corresponding $y$ value and area are found by:
```{julia}
xstar = soln[]
height(xstar), A(xstar)
```
The `minimum` property also holds the identified minimum:
```{julia}
soln.minimum # compare with A(soln[], nothing)
```
The package is a wrapper around other packages. The output of the underlying package is presented in the `original` property:
```{julia}
soln.original
```
---
This problem can also be approached symbolically, using `ModelingToolkit`.
For example, we set up the problem with:
```{julia}
using ModelingToolkit
@parameters P
@variables x
y = (P - 2x)/2
Area = - x*y
```
The above should be self explanatory. To put into a form to pass to `solve` we define a "system" by specifying our objective function, the variables, and the parameters.
```{julia}
@named sys = OptimizationSystem(Area, [x], [P])
```
(This step is different, as before an `OptimizationFunction` was defined; we use `@named`, as above, to ensure the system has the same name as the identifier, `sys`.)
This system is passed to `OptimizationProblem` along with a specification of the initial condition ($x=4$) and the perimeter ($P=25$). A vector of pairs is used below:
```{julia}
prob = OptimizationProblem(sys, [x => 4.0], [P => 25.0]; grad=true, hess=true)
```
The keywords `grad=true` and `hess=true` instruct for automatic derivatives to be taken as needed. These are needed in the choice of method, `Newton`, below.
Solving this problem then follows the same pattern as before, again with `Newton` we have:
```{julia}
solve(prob, Newton())
```
(A derivative-free method like `NelderMead()` could be used and then the `grad` and `hess` keywords above would be unnecessary, though not harmful.)
---
The related calculus problem:
> Among all rectangles with a fixed area, find the one with *minimum* perimeter
could be similarly approached:
```{julia}
@parameters Area
@variables x
y = Area/x # from A = xy
P = 2x + 2y
@named sys = OptimizationSystem(P, [x], [Area])
u0 = [x => 4.0]
p = [Area => 25.0]
prob = OptimizationProblem(sys, u0, p; grad=true, hess=true)
soln = solve(prob, LBFGS())
```
We used an initial guess of $x=4$ above. The `LBFGS` method is a computationally efficient modification of the Broyden-Fletcher-Goldfarb-Shanno algorithm ... It is a quasi-Newton method that updates an approximation to the Hessian using past approximations as well as the gradient." On this problem it performs similarly to `Newton`, though in general may be preferable for higher-dimensional problems.
### Two dimensional
Scalar functions of two input variables can have their minimum value identified in the same manner using `Optimization.jl`.
For example, consider the function
$$
f(x,y) = (x + 2y - 7)^2 + (2x + y - 5)^2
$$
We wish to minimize this function.
We begin by defining a function in `Julia`:
```{julia}
function f(u, p)
x, y = u
(x + 2y - 7)^2 + (2x + y - 5)^2
end
```
We turn this into an optimization function by specifying how derivatives will be taken, as we will the `LBFGS` algorithm below:
```{julia}
ff = OptimizationFunction(f, Optimization.AutoForwardDiff())
```
We will begin our search at the origin. We have to be mindful to use floating point numbers here:
```{julia}
u0 = [0.0, 0.0] # or zeros(2)
```
```{julia}
prob = OptimizationProblem(ff, u0)
```
Finally, we solve the values:
```{julia}
solve(prob, LBFGS())
```
The value of $(1, 3)$ agrees with the contour graph, as it is a point in the interior of the contour for the smallest values displayed.
```{julia}
using Plots
xs = range(0, 2, length=100)
ys = range(2, 4, length=100)
contour(xs, ys, (x,y) -> f((x,y), nothing))
```
We could also use a *derivative-free* method, and skip a step:
```{julia}
prob = OptimizationProblem(f, u0) # skip making an OptimizationFunction
solve(prob, NelderMead())
```
## Integration (`Integrals.jl`)
The `Integrals` package provides a common interface to different numeric integration packages in the `Julia` ecosystem. For example, `QuadGK` and `HCubature`. The value of this interface, over those two packages, is its non-differentiated access to other packages, which for some uses may be more performant.
The package follows the same `problem-algorithm-solve` interface, as already seen.
The interface is designed for $1$-and-higher dimensional integrals.
The package is loaded with
```{julia}
using Integrals
```
For a simple definite integral, such as $\int_0^\pi \sin(x)dx$, we have:
```{julia}
f(x, p) = sin(x)
prob = IntegralProblem(f, 0.0, pi)
soln = solve(prob, QuadGKJL())
```
To get access to the answer, we can use indexing notation:
```{julia}
soln[]
```
Comparing to just using `QuadGK`, the same definite integral would be estimated with:
```{julia}
using QuadGK
quadgk(sin, 0, pi)
```
The estimated upper bound on the error from `QuadGK`, is available through the `resid` property on the `Integrals` output:
```{julia}
soln.resid
```
The `Integrals` solution is a bit more verbose, but it is more flexible. For example, the `HCubature` package provides a similar means to compute $n$- dimensional integrals. For this problem, the modifications would be:
```{julia}
f(x, p) = sin.(x)
prob = IntegralProblem(f, [0.0], [pi])
soln = solve(prob, HCubatureJL())
```
```{julia}
soln[]
```
The estimated maximum error is also given by `resid`:
```{julia}
soln.resid
```
---
As well, suppose we wanted to parameterize our function and then differentiate.
Consider $d/dp \int_0^\pi \sin(px) dx$. We can do this integral directly to get
$$
\begin{align*}
\frac{d}{dp} \int_0^\pi \sin(px)dx
&= \frac{d}{dp}\left( \frac{-1}{p} \cos(px)\Big\rvert_0^\pi\right)\\
&= \frac{d}{dp}\left( -\frac{\cos(p\cdot\pi)-1}{p}\right)\\
&= \frac{\cos(p\cdot \pi) - 1)}{p^2} + \frac{\pi\cdot\sin(p\cdot\pi)}{p}
\end{align*}
$$
Using `Integrals` with `QuadGK` we have:
```{julia}
f(x, p) = sin(p*x)
function ∫sinpx(p)
prob = IntegralProblem(f, 0.0, pi, p)
solve(prob, QuadGKJL())
end
```
We can compute values at both $p=1$ and $p=2$:
```{julia}
∫sinpx(1), ∫sinpx(2)
```
To find the derivative in $p$ , we have:
```{julia}
ForwardDiff.derivative(∫sinpx, 1), ForwardDiff.derivative(∫sinpx, 2)
```
(In `QuadGK`, the following can be differentiated `∫sinpx(p) = quadgk(x -> sin(p*x), 0, pi)[1]` as well. `Integrals` gives a consistent interface.
### Higher dimension integrals
The power of a common interface is the ability to swap backends and the uniformity for different dimensions. Here we discuss integrals of scalar-valued and vector-valued functions.
#### $f: R^n \rightarrow R$
The area under a surface generated by $z=f(x,y)$ over a rectangular region $[a,b]\times[c,d]$ can be readily computed. The two coding implementations require $f$ to be expressed as a function of a vector–*and* a parameter–and the interval to be expressed using two vectors, one for the left endpoints (`[a,c]`) and on for the right endpoints (`[b,d]`).
For example, the area under the function $f(x,y) = 1 + x^2 + 2y^2$ over $[-1/2, 1/2] \times [-1,1]$ is computed by:
```{julia}
f(x, y) = 1 + x^2 + 2y^2 # match math
fxp(x, p) = f(x[1], x[2]) # prepare for IntegralProblem
ls = [-1/2, -1] # left endpoints
rs = [1/2, 1] # right endpoints
prob = IntegralProblem(fxp, ls, rs)
soln = solve(prob, HCubatureJL())
```
Of course, we could have directly defined the function (`fxp`) using indexing of the `x` variable.
---
For non-rectangular domains a change of variable is required.
For example, an integral to assist in finding the volume of a sphere might be
$$
V = 2 \iint_R \sqrt{\rho^2 - x^2 - y^2} dx dy
$$
where $R$ is the disc of radius $\rho$ in the $x-y$ plane.
The usual approach is to change to polar-coordinates and write this integral as
$$
V = \int_0^{2\pi}\int_0^\rho \sqrt{\rho^2 - r^2} r dr d\theta
$$
the latter being an integral over a rectangular domain.
To compute this transformed integral, we might have:
```{julia}
function vol_sphere(ρ)
f(rθ, p) = sqrt(ρ^2 - rθ[1]^2) * rθ[1]
ls = [0,0]
rs = [ρ, 2pi]
prob = IntegralProblem(f, ls, rs)
solve(prob, HCubatureJL())
end
vol_sphere(2)
```
If it is possible to express the region to integrate as $G(R)$ where $R$ is a rectangular region, then the change of variables formula,
$$
\iint_{G(R)} f(x) dA = \iint_R (f\circ G)(u) |det(J_G(u)| dU
$$
turns the integral into the non-rectangular domain $G(R)$ into one over the rectangular domain $R$. The key is to *identify* $G$ and to compute the Jacobian. The latter is simply accomplished with `ForwardDiff.jacobian`.
For an example, we find the moment of inertia about the axis of the unit square tilted counter-clockwise an angle $0 \leq \alpha \leq \pi/2$.
The counter clockwise rotation of a unit square by angle $\alpha$ is described by
$$
G(u, v) = \langle \cos(\alpha)\cdot u - \sin(\alpha)\cdot v, \sin(\alpha)\cdot u, +\cos(\alpha)\cdot v \rangle
$$
So we have $\iint_{G(R)} x^2 dA$ is computed by the following with $\alpha=\pi/4$:
```{julia}
import LinearAlgebra: det
𝑓(uv) = uv[1]^2
function G(uv)
α = pi/4 # could be made a parameter
u,v = uv
[cos(α)*u - sin(α)*v, sin(α)*u + cos(α)*v]
end
f(u, p) = (𝑓∘G)(u) * det(ForwardDiff.jacobian(G, u))
prob = IntegralProblem(f, [0,0], [1,1])
solve(prob, HCubatureJL())
```
#### $f: R^n \rightarrow R^m$
The `Integrals` package provides an interface for vector-valued functions. By default, the number of dimensions in the output is assumed to be $1$, but the `nout` argument can adjust that.
Let $f$ be vector valued with components $f_1, f_2, \dots, f_m$, then the output below is the vector with components $\iint_R f_1 dV, \iint_R f_2 dV, \dots, \iint_R f_m dV$.
For a trivial example, we have:
```{julia}
f(x, p) = [x[1], x[2]^2]
prob = IntegralProblem(f, [0,0],[3,4], nout=2)
solve(prob, HCubatureJL())
```

10
quarto/alternatives/sciml_interface.__jmd__
View File

@@ -1,10 +0,0 @@
# seems like with
* quadrature (renamed)
* nonlinsolve
* GalacticOptim (renamed)
* symbolic-numeric integration
* symbolics.jl
...
This should be mentioned

392
quarto/alternatives/symbolics.__jmd__
View File

@@ -1,392 +0,0 @@
# Symbolics.jl
Incorporate:
Basics
https://github.com/SciML/ModelingToolkit.jl
https://github.com/JuliaSymbolics/Symbolics.jl
https://github.com/JuliaSymbolics/SymbolicUtils.jl
* Rewriting
https://github.com/JuliaSymbolics/SymbolicUtils.jl
* Plotting
Polynomials
Limits
XXX ... room here!
Derivatives
https://github.com/JuliaSymbolics/Symbolics.jl
Integration
https://github.com/SciML/SymbolicNumericIntegration.jl
The `Symbolics.jl` package is a Computer Algebra System (CAS) built entirely in `Julia`.
This package is under heavy development.
## Algebraic manipulations
### construction
@variables
SymbolicUtils.@syms assumptions
x is a `Num`, `Symbolics.value(x)` is of type `SymbolicUtils{Real, Nothing}
relation to SymbolicUtils
Num wraps things; Term
### Substitute
### Simplify
simplify
expand
rewrite rules
### Solving equations
solve_for
## Expressions to functions
build_function
## Derivatives
1->1: Symbolics.derivative(x^2 + cos(x), x)
1->3: Symbolics.derivative.([x^2, x, cos(x)], x)
3 -> 1: Symbolics.gradient(x*y^z, [x,y,z])
2 -> 2: Symbolics.jacobian([x,y^z], [x,y])
# higher order
1 -> 1: D(ex, x, n=1) = foldl((ex,_) -> Symbolics.derivative(ex, x), 1:n, init=ex)
2 -> 1: (2nd) Hessian
## Differential equations
## Integrals
WIP
## ----
# follow sympy tutorial
using Symbolics
import SymbolicUtils
@variables x y z
# substitution
ex = cos(x) + 1
substitute(ex, Dict(x=>y))
substitute(ex, Dict(x=>0)) # does eval
ex = x^y
substitute(ex, Dict(y=> x^y))
# expand trig
r1 = @rule sin(2 * ~x) => 2sin(~x)*cos(~x)
r2 = @rule cos(2 * ~x) => cos(~x)^2 - sin(~x)^2
expand_trig(ex) = simplify(ex, RuleSet([r1, r2]))
ex = sin(2x) + cos(2x)
expand_trig(ex)
## Multiple
@variables x y z
ex = x^3 + 4x*y -z
substitute(ex, Dict(x=>2, y=>4, z=>0))
# Converting Strings to Expressions
# what is sympify?
# evalf
# lambdify: symbolic expression -> function
ex = x^3 + 4x*y -z
λ = build_function(ex, x,y,z, expression=Val(false))
λ(2,4,0)
# pretty printing
using Latexify
latexify(ex)
# Simplify
@variables x y z t
simplify(sin(x)^2 + cos(x)^2)
simplify((x^3 + x^2 - x - 1) / (x^2 + 2x + 1)) # fails, no factor
simplify(((x+1)*(x^2-1))/((x+1)^2)) # works
import SpecialFunctions: gamma
simplify(gamma(x) / gamma(x-2)) # fails
# Polynomial
## expand
expand((x+1)^2)
expand((x+2)*(x-3))
expand((x+1)*(x-2) - (x-1)*x)
## factor
### not defined
## collect
COLLECT_RULES = [
@rule(~x*x^(~n::SymbolicUtils.isnonnegint) => (~x, ~n))
@rule(~x * x => (~x, 1))
]
function _collect(ex, x)
d = Dict()
exs = expand(ex)
if SymbolicUtils.operation(Symbolics.value(ex)) != +
d[0] => ex
else
for aᵢ ∈ SymbolicUtils.arguments(Symbolics.value(expand(ex)))
u = simplify(aᵢ, RuleSet(COLLECT_RULES))
if isa(u, Tuple)
a,n = u
else
a,n = u,0
end
d[n] = get(d, n, 0) + a
end
end
d
end
## cancel -- no factor
## apart -- no factor
## Trignometric simplification
INVERSE_TRIG_RUELS = [@rule(cos(acos(~x)) => ~x)
@rule(acos(cos(~x)) => abs(rem2pi(~x, RoundNearest)))
@rule(sin(asin(~x)) => ~x)
@rule(asin(sin(~x)) => abs(rem2pi(x + pi/2, RoundNearest)) - pi/2)
]
@variables θ
simplify(cos(acos(θ)), RuleSet(INVERSE_TRIG_RUELS))
# Copy from https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/master/src/simplify_rules.jl
# the TRIG_RULES are applied by simplify by default
HTRIG_RULES = [
@acrule(-sinh(~x)^2 + cosh(~x)^2 => one(~x))
@acrule(sinh(~x)^2 + 1 => cosh(~x)^2)
@acrule(cosh(~x)^2 + -1 => -sinh(~x)^2)
@acrule(tanh(~x)^2 + 1*sech(~x)^2 => one(~x))
@acrule(-tanh(~x)^2 + 1 => sech(~x)^2)
@acrule(sech(~x)^2 + -1 => -tanh(~x)^2)
@acrule(coth(~x)^2 + -1*csch(~x)^2 => one(~x))
@acrule(coth(~x)^2 + -1 => csch(~x)^2)
@acrule(csch(~x)^2 + 1 => coth(~x)^2)
@acrule(tanh(~x) => sinh(~x)/cosh(~x))
@acrule(sinh(-~x) => -sinh(~x))
@acrule(cosh(-~x) => -cosh(~x))
]
trigsimp(ex) = simplify(simplify(ex, RuleSet(HTRIG_RULES)))
trigsimp(sin(x)^2 + cos(x)^2)
trigsimp(sin(x)^4 -2cos(x)^2*sin(x)^2 + cos(x)^4) # no factor
trigsimp(cosh(x)^2 + sinh(x)^2)
trigsimp(sinh(x)/tanh(x))
EXPAND_TRIG_RULES = [
@acrule(sin(~x+~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y))
@acrule(sinh(~x+~y) => sinh(~x)*cosh(~y) + cosh(~x)*sinh(~y))
@acrule(sin(2*~x) => 2sin(~x)*cos(~x))
@acrule(sinh(2*~x) => 2sinh(~x)*cosh(~x))
@acrule(cos(~x+~y) => cos(~x)*cos(~y) - sin(~x)*sin(~y))
@acrule(cosh(~x+~y) => cosh(~x)*cosh(~y) + sinh(~x)*sinh(~y))
@acrule(cos(2*~x) => cos(~x)^2 - sin(~x)^2)
@acrule(cosh(2*~x) => cosh(~x)^2 + sinh(~x)^2)
@acrule(tan(~x+~y) => (tan(~x) - tan(~y)) / (1 + tan(~x)*tan(~y)))
@acrule(tanh(~x+~y) => (tanh(~x) + tanh(~y)) / (1 + tanh(~x)*tanh(~y)))
@acrule(tan(2*~x) => 2*tan(~x)/(1 - tan(~x)^2))
@acrule(tanh(2*~x) => 2*tanh(~x)/(1 + tanh(~x)^2))
]
expandtrig(ex) = simplify(simplify(ex, RuleSet(EXPAND_TRIG_RULES)))
expandtrig(sin(x+y))
expandtrig(tan(2x))
# powers
# in genearl x^a*x^b = x^(a+b)
@variables x y a b
simplify(x^a*x^b - x^(a+b)) # 0
# x^a*y^a = (xy)^a When x,y >=0, a ∈ R
simplify(x^a*y^a - (x*y)^a)
## ??? How to specify such assumptions?
# (x^a)^b = x^(ab) only if b ∈ Int
@syms x a b
simplify((x^a)^b - x^(a*b))
@syms x a b::Int
simplify((x^a)^b - x^(a*b)) # nope
ispositive(x) = isa(x, Real) && x > 0
_isinteger(x) = isa(x, Integer)
_isinteger(x::SymbolicUtils.Sym{T,S}) where {T <: Integer, S} = true
POWSIMP_RULES = [
@acrule((~x::ispositive)^(~a::isreal) * (~y::ispositive)^(~a::isreal) => (~x*~y)^~a)
@rule(((~x)^(~a))^(~b::_isinteger) => ~x^(~a * ~b))
]
powsimp(ex) = simplify(simplify(ex, RuleSet(POWSIMP_RULES)))
@syms x a b::Int
simplify((x^a)^b - x^(a*b)) # nope
EXPAND_POWER_RULES = [
@rule((~x)^(~a + ~b) => (_~)^(~a) * (~x)^(~b))
@rule((~x*~y)^(~a) => (~x)^(~a) * (~y)^(~a))
## ... more on simplification...
## Calculus
@variables x y z
import Symbolics: derivative
derivative(cos(x), x)
derivative(exp(x^2), x)
# multiple derivative
Symbolics.derivative(ex, x, n::Int) = reduce((ex,_) -> derivative(ex, x), 1:n, init=ex) # helper
derivative(x^4, x, 3)
ex = exp(x*y*z)
using Chain
@chain ex begin
derivative(x, 3)
derivative(y, 3)
derivative(z, 3)
end
# using Differential operator
expr = exp(x*y*z)
expr |> Differential(x)^2 |> Differential(y)^3 |> expand_derivatives
# no integrate
# no limit
# Series
function series(ex, x, x0=0, n=5)
Σ = zero(ex)
for i ∈ 0:n
ex = expand_derivatives((Differential(x))(ex))
Σ += substitute(ex, Dict(x=>0)) * x^i / factorial(i)
end
Σ
end
# finite differences
# Solvers
@variables x y z a
eq = x ~ a
Symbolics.solve_for(eq, x)
eqs = [x + y + z ~ 1
x + y + 2z ~ 3
x + 2y + 3z ~ 3
]
vars = [x,y,z]
xs = Symbolics.solve_for(eqs, vars)
[reduce((ex, r)->substitute(ex, r), Pair.(vars, xs), init=ex.lhs) for ex ∈ eqs] == [eq.rhs for eq ∈ eqs]
A = [1 1; 1 2]
b = [1, 3]
xs = Symbolics.solve_for(A*[x,y] .~ b, [x,y])
A*xs - b
A = [1 1 1; 1 1 2]
b = [1,3]
A*[x,y,z] - b
Symbolics.solve_for(A*[x,y,z] .~ b, [x,y,z]) # fails, singular
# nonlinear solve
# use `λ = mk_function(ex, args, expression=Val(false))`
# polynomial roots
# differential equations

201
quarto/alternatives/symbolics.qmd
View File

@@ -9,7 +9,7 @@ There are a few options in `Julia` for symbolic math, for example, the `SymPy` p
## About
The `Symbolics` package bill itself as a "fast and modern Computer Algebra System (CAS) for a fast and modern programming language." This package relies on the `SymbolicUtils` package and is built upon by the `ModelingToolkit` package, which we don't describe here.
The `Symbolics` package bills itself as a "fast and modern Computer Algebra System (CAS) for a fast and modern programming language." This package relies on the `SymbolicUtils` package and is built upon by the `ModelingToolkit` package, which we don't describe here.
We begin by loading the `Symbolics` package which when loaded re-exports the `SymbolicUtils` package.
@@ -22,7 +22,7 @@ using Symbolics
## Symbolic variables
Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two package discussed here.
Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two packages discussed here.
`SymbolicUtils` creates variables which carry `Julia` type information (e.g. `Int`, `Float64`, ...). This type information carries through operations involving these variables. Symbolic variables can be created with the `@syms` macro. For example
@@ -32,10 +32,10 @@ Symbolic math at its core involves symbolic variables, which essentially defer e
@syms x y::Int f(x::Real)::Real
```
This creates `x` a symbolic value with type `Number`, `y` a symbolic variable holding integer values, and `f` a symbolic function of a single real variable outputing a real variable.
This creates `x` a symbolic value with symbolic type `Number`, `y` a symbolic variable holding integer values, and `f` a symbolic function of a single real variable outputting a real variable.
The `symtype` function reveals the underlying type:
The non-exported `symtype` function reveals the underlying type:
```{julia}
@@ -44,14 +44,14 @@ import Symbolics.SymbolicUtils: symtype
symtype(x), symtype(y)
```
For `y`, the symbolic type being real does not imply the `y` has a subtype of `Real`:
For `y`, the symbolic type being real does not imply the type of `y` is a subtype of `Real`:
```{julia}
isa(y, Real)
```
We see that the function `f`, when called with `y` would return a value of (symbolic) type `Real`:
We see that the function `f` when called with `y` would return a value of (symbolic) type `Real`:
```{julia}
@@ -68,7 +68,7 @@ f(x)
:::{.callout-note}
## Note
The `SymPy` package also has an `@syms` macro to create variables. Though they names agree, they do different things. Using both packages together would require qualifying many shared method names. For `SymbolicUtils`, the `@syms` macro uses `Julia` types to parameterize the variables. In `SymPy` it is possible to specify *assumptions* on the variables, but that is different and not useful for dispatch without some extra effort.
The `SymPy` package also has an `@syms` macro to create variables. Though their names agree, they do different things. Using both packages together would require qualifying many shared method names. For `SymbolicUtils`, the `@syms` macro uses `Julia` types to parameterize the variables. In `SymPy` it is possible to specify *assumptions* on the variables, but that is different and not useful for dispatch without some extra effort.
:::
@@ -137,7 +137,7 @@ typeof(sin(x)), typeof(Symbolics.value(sin(x)))
### Tree structure to expressions
The `TermInterface` package is used by `SymbolicUtils` to explore the tree structdure of an expression. The main methods are (cf. [symbolicutils.jl](https://symbolicutils.juliasymbolics.org/#expression_interface)):
The `TermInterface` package is used by `SymbolicUtils` to explore the tree structure of an expression. The main methods are (cf. [SymbolicUtils.jl](https://symbolicutils.juliasymbolics.org/#expression_interface)):
* `istree(ex)`: `true` if `ex` is not a *leaf* node (like a symbol or numeric literal)
@@ -221,7 +221,7 @@ free_symbols(ex)
### Substitute
The `substitute` command is used to replace values with other values. For examples:
The `substitute` command is used to replace values with other values. For example:
```{julia}
@@ -260,7 +260,7 @@ substitute(ex, x=>π), substitute(ex, x=>π, fold=false)
Algebraic operations with symbolic values can involve an exponentially increasing number of terms. As such, some simplification rules are applied after an operation to reduce the complexity of the computed value.
For example, `0+x` should simplify to `x`, `0+x`, `x^0`, or `x^1` should simplify, to some natural answer.
For example, `0+x` should simplify to `x`, as well `1*x`, `x^0`, or `x^1` should each simplify, to some natural answer.
`SymbolicUtils` also [simplifies](https://symbolicutils.juliasymbolics.org/#simplification) several other expressions, including:
@@ -269,7 +269,7 @@ For example, `0+x` should simplify to `x`, `0+x`, `x^0`, or `x^1` should simplif
* `-x` becomes `(-1)*x`
* `x * x` becomes `x^2` (and `x^n` if more terms). Meaning this expression is represented as a power, not a product
* `x + x` becomes `2*x` (and `n*x` if more terms). Similarly, this represented as a product, not a sum.
* `p/q * x` becomes `(p*x)/q)`, similarly `p/q * x/y` becomes `(p*x)/(q*y)`
* `p/q * x` becomes `(p*x)/q)`, similarly `p/q * x/y` becomes `(p*x)/(q*y)`. (Division wraps multiplication.)
In `SymbolicUtils`, this *rewriting* is accomplished by means of *rewrite rules*. The package makes it easy to apply user-written rewrite rules.
@@ -278,7 +278,7 @@ In `SymbolicUtils`, this *rewriting* is accomplished by means of *rewrite rules*
### Rewriting
Many algebraic simplifications are done by the `simplify` command. For example, the basic trignometric identities are applied:
Many algebraic simplifications are done by the `simplify` command. For example, the basic trigonometric identities are applied:
```{julia}
@@ -307,11 +307,14 @@ ex1 = substitute(ex, x => sin(x + y + z))
ex1 |> Symbolics.value |> r |> Num
```
Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trignometry and exponential simplifications:
Rewrite rules when applied return the rewritten expression, if there is a match, or `nothing`.
Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trigonometry and exponential simplifications:
```{julia}
r = @rule(exp(~x)^(~y) => exp(~x * ~y))
r = @rule(exp(~x)^(~y) => exp(~x * ~y)) # (e^x)^y -> e^(x*y)
ex = exp(-x+z)^y
ex, ex |> Symbolics.value |> r |> Num
```
@@ -346,7 +349,7 @@ ex = exp(-x + z)^y
Symbolics.toexpr(ex)
```
This output shows an internal representation of the steps for computing the value `ex` given different inputs.
This output shows an internal representation of the steps for computing the value `ex` given different inputs. (The number `(-1)` multiplies `x`, this is added to `z` and the result passed to `exp`. That values is then used as the base for `^` with exponent `y`.)
Such `Julia` expressions are one step away from building `Julia` functions for evaluating symbolic expressions fast (though with some technical details about "world age" to be reckoned with). The `build_function` function with the argument `expression=Val(false)` will compile a `Julia` function:
@@ -376,7 +379,7 @@ The documentation colorfully says "`build_function` is kind of like if `lambdify
The above, through passing $3$ variables after the expression, creates a function of $3$ variables. Functions of a vector of inputs can also be created, just by expressing the variables in that manner:
```{juila}
```{julia}
h1 = build_function(ex, [x, y, z]; expression=Val(false))
h1([1, 2, 3]) # not h1(1,2,3)
```
@@ -398,7 +401,7 @@ Roots.find_zero(λ, (1, 2))
### Plotting
Using `Plots`, the plotting symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via `build_function`.
Using `Plots`, the plotting of symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via `build_function`.
For example,
@@ -451,7 +454,7 @@ For example
d, r = polynomial_coeffs(a*x^2 + b*x + c, (x,))
```
The first term output is dictionary who's keys are the monomials and who's values are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just the constant `c`.
The first term output is dictionary with keys which are the monomials and with values which are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just the constant `c`.
The expression can then be reconstructed through
@@ -461,7 +464,7 @@ The expression can then be reconstructed through
r + sum(v*k for (k,v) ∈ d)
```
The above has `a,b,c` as parameters and `x` as the symbol. This separation is designated by passing the desired polynomials symbols to `polynomial_coeff` as an iterable. (Above as a $1$-element tuple.)
The above has `a,b,c` as parameters and `x` as the symbol. This separation is designated by passing the desired polynomial symbols to `polynomial_coeff` as an iterable. (Above as a $1$-element tuple.)
More complicated polynomials can be similarly decomposed:
@@ -552,6 +555,69 @@ m,n = degree.(nd(ex))
m > n ? "limit is infinite" : m < n ? "limit is 0" : "limit is a constant"
```
### Vectors and matrices
Symbolic vectors and matrices can be created with a specified size:
```{julia}
@variables v[1:3] M[1:2, 1:3] N[1:3, 1:3]
```
Computations, like finding the determinant below, are lazy unless the values are `collect`ed:
```{julia}
using LinearAlgebra
det(N)
```
```{julia}
det(collect(N))
```
Similarly, with `norm`:
```{julia}
norm(v)
```
and
```{julia}
norm(collect(v))
```
Matrix multiplication is also deferred, but the size compatability of the matrices and vectors is considered early:
```{julia}
M*N, N*N, M*v
```
This errors, as the matrix dimensions are not compatible for multiplication:
```{julia}
#| error: true
N*M
```
Similarly, linear solutions can be symbolically specified:
```{julia}
@variables R[1:2, 1:2] b[1:2]
R \ b
```
```{julia}
collect(R \ b)
```
### Algebraically solving equations
@@ -567,14 +633,33 @@ or
```{julia}
ex = [5x + 2y, 6x + 3y] .~ [1, 2]
eqs = [5x + 2y, 6x + 3y] .~ [1, 2]
```
The `Symbolics.solve_for` function can solve *linear* equations. For example,
```{julia}
Symbolics.solve_for(ex, [x, y])
Symbolics.solve_for(eqs, [x, y])
```
The coefficients can be symbolic. Two examples could be:
```{julia}
@variables m b x y
eq = y ~ m*x + b
Symbolics.solve_for(eq, x)
```
```{julia}
@variables a11 a12 a22 x y b1 b2
R,X,b = [a11 a12; 0 a22], [x; y], [b1, b2]
eqs = R*X .~ b
```
```{julia}
Symbolics.solve_for(eqs, [x,y])
```
### Limits
@@ -591,26 +676,43 @@ As of writing, there is no extra functionality provided by `Symbolics` for compu
```{julia}
@variables a b c x
ex = a*x^2 + b*x + c
Symbolics.derivative(ex, x)
y = a*x^2 + b*x + c
yp = Symbolics.derivative(y, x)
```
The computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:
Or to find a critical point:
```{julia}
Symbolics.solve_for(yp ~ 0, x) # linear equation to solve
```
The derivative computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:
```{julia}
D = Differential(x)
D(ex)
D(y)
```
and then
```{julia}
expand_derivatives(D(ex))
expand_derivatives(D(y))
```
The differentials can be multiplied to create operators for taking higher-order derivatives:
Using `Differential`, differential equations can be specified. An example was given in [ODEs](../ODEs/differential_equations.html), using `ModelingToolkit`.
Higher order derivatives can be done through composition:
```{julia}
D(D(y)) |> expand_derivatives
```
Differentials can also be multiplied to create operators for taking higher-order derivatives:
```{julia}
@@ -628,7 +730,7 @@ In addition to `Symbolics.derivative` there are also the helper functions, such
Symbolics.hessian(ex, [x,y])
```
The `gradient` function is also available
The `gradient` function is also defined
```{julia}
@@ -637,13 +739,13 @@ ex = x^2 - 2x*y + z*y
Symbolics.gradient(ex, [x, y, z])
```
The `jacobian` takes an array of expressions:
The `jacobian` function takes an array of expressions:
```{julia}
@variables x y
exs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(exs, [x,y])
eqs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(eqs, [x,y])
```
### Integration
@@ -655,13 +757,15 @@ The `SymbolicNumericIntegration` package provides a means to integrate *univaria
Symbolic integration can be approached in different ways. SymPy implements part of the Risch algorithm in addition to other algorithms. Rules-based algorithms could also be implemented.
For example, here is a simple rule that could be used to integrate a single integral
For a trivial example, here is a rule that could be used to integrate a single integral
```{julia}
is_var(x) = (xs = Symbolics.get_variables(x); length(xs) == 1 && xs[1] === x)
@syms x ∫(x)
is_var(x) = (xs = Symbolics.get_variables(x); length(xs) == 1 && xs[1] === x)
r = @rule ∫(~x::is_var) => x^2/2
r(∫(x))
```
@@ -671,23 +775,37 @@ The `SymbolicNumericIntegration` package includes many more predicates for doing
If $f(x)$ is to be integrated, a set of *candidate* answers is generated. The following is **proposed** as an answer: $\sum q_i \Theta_i(x)$. Differentiating the proposed answer leads to a *linear system of equations* that can be solved.
The example in the [paper](https://arxiv.org/pdf/2201.12468v2.pdf) describing the method is with $f(x) = x \sin(x)$ and the candidate thetas are ${x, \sin(x), \cos(x), x\sin(x), x\cos(x)}$ so that we propose:
The example in the [paper](https://arxiv.org/pdf/2201.12468v2.pdf) describing the method is with $f(x) = x \sin(x)$ and the candidate thetas are ${x, \sin(x), \cos(x), x\sin(x), x\cos(x)}$ so that the propose answer is:
$$
\int f(x) dx = q_1 x + q_2 \sin(x) + q_3 \cos(x) + q_4 x \sin(x) + q_4 x \cos(x)
$$
Differentiating both sides, yields a term $x\sin(x)$ on the left, and equating coefficients gives:
We differentiate the right hand side:
```{julia}
@variables q[1:5] x
ΣqᵢΘᵢ = dot(collect(q), (x, sin(x), cos(x), x*sin(x), x*cos(x)))
simplify(Symbolics.derivative(ΣqᵢΘᵢ, x))
```
This must match $x\sin(x)$ so we have by equating coefficients of the respective terms:
$$
q_1 = q_4 = 0,\quad q_5 = -1, \quad q_4 - q_3 = q_2 - q_5 = 0
q_2 + q_5 = 0, \quad q_4 = 0, \quad q_1 = 0, \quad q_3 = 0, \quad q_5 = -1
$$
which can be solved with $q_5=-1$, $q_2=1$, and the other coefficients being $0$. That is $\int f(x) dx = 1 \sin(x) + (-1) x\cos(x)$.
That is $q_2=1$, $q_5=-1$, and the other coefficients are $0$, giving an answer computed with:
```{julia}
d = Dict(q[i] => v for (i,v) ∈ enumerate((0,1,0,0,-1)))
substitute(ΣqᵢΘᵢ, d)
```
The package provides an algorithm for the creation of candidates and the means to solve when possible. The `integrate` function is the main entry point. It returns three values: `solved`, `unsolved`, and `err`. The `unsolved` is the part of the integrand which can not be solved through this package. It is `0` for a given problem when `integrate` is successful in identifying an antiderivative, in which case `solved` is the answer. The value of `err` is a bound on the numerical error introduced by the algorithm.
@@ -701,7 +819,7 @@ using SymbolicNumericIntegration
integrate(x * sin(x))
```
The second term is `0`, as this has an identified antiderivative.
The second term is `0`, as this integrand has an identified antiderivative.
```{julia}
@@ -725,6 +843,9 @@ The derivative of `u` matches up to some numeric tolerance:
Symbolics.derivative(u, x) - sin(x)^5
```
---
The integration of rational functions (ratios of polynomials) can be done algorithmically, provided the underlying factorizations can be identified. The `SymbolicNumericIntegration` package has a function `factor_rational` that can identify factorizations.
@@ -758,14 +879,14 @@ u = 1 / expand((x^2+1)*(x-2)^2)
v = factor_rational(u)
```
As such, the integrals have numeric differences:
As such, the integrals have numeric differences from their mathematical counterparts:
```{julia}
a,b,c = integrate(u)
```
We can see a bit of why through the following which needs a tolerance set to identify the rational numbers correctly:
We can see a bit of how much through the following, which needs a tolerance set to identify the rational numbers of the mathematical factorization correctly:
```{julia}

0
.github/workflows/documentation.yml → quarto/documentation.yml
View File

2
quarto/index.qmd
View File

@@ -74,7 +74,7 @@ being allowed per notebook, there is frequent use of
[Unicode](./misc/unicode.html) symbols for variable names.
To contribute -- say by suggesting addition topics, correcting a
mistake, or fixing a typo -- click the "Edit this page" link.
mistake, or fixing a typo -- click the "Edit this page" link and join the list of [contributors](https://github.com/jverzani/CalculusWithJuliaNotes.jl/graphs/contributors).
----

8
quarto/precalc/functions.qmd
View File

@@ -199,7 +199,7 @@ end
The conditions for the `if` statements are expressions that evaluate to either `true` or `false`, such as generated by the Boolean operators `<`, `<=`, `==`, `!-`, `>=`, and `>`.
If familiar with `if` conditions, they are natural to use. However, for simpler cases of "if-else" `Julia` provides the more convenient *ternary* operator: `cond ? if_true : if_false`. (The name comes from the fact that there are three arguments specified.) The ternary operator checks the condition and if true returns the first expression, whereas if the condition is false the second condition is returned. Both expressions are evaluated. (The [short-circuit](https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation) operators can be used to avoid both evaluations.)
If familiar with `if` conditions, they are natural to use. However, for simpler cases of "if-else" `Julia` provides the more convenient *ternary* operator: `cond ? if_true : if_false`. (The name comes from the fact that there are three arguments specified.) The ternary operator checks the condition and if true returns the first expression, whereas if the condition is false the second condition is returned. (Another useful control flow construct is [short-circuit](https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation) evaluation.)
For example, here is one way to define an absolute value function:
@@ -1129,6 +1129,11 @@ answ = 2
radioq(choices, answ, keep_order=true)
```
:::{.callout-note}
## Note
The parentheses in `(sin ∘ cos)(pi/4)` are needed due to the order of operations, with `cos(pi/4)` being evaluated first in the expression `sin ∘ cos(pi/4)`. Alternatively, one can define a function `sc = sin ∘ cos` (without parentheses), then call it through `sc(pi/4)`.
:::
###### Question
@@ -1344,4 +1349,3 @@ The interval is a nearly exact estimate, as guaranteed by `IntervalArithmetic`.
"""]
radioq(choices, 1)
```

4
quarto/precalc/julia_overview.qmd
View File

@@ -321,7 +321,7 @@ log
### Built-in functions
`Julia` has numerous built-in [mathematical](http://julia.readthedocs.io/) functions, we review a few here:
`Julia` has numerous built-in [mathematical](https://docs.julialang.org/en/v1/manual/mathematical-operations/) functions, we review a few here:
#### Powers logs and roots
@@ -631,5 +631,3 @@ plot(64 - (1/2)*32 * x^2, 0, 2)
* SymPy has functions for manipulating expressions: `simplify`, `expand`, `together`, `factor`, `cancel`, `apart`, $...$
* SymPy has functions for basic math: `factor`, `roots`, `solve`, `solveset`, $\dots$
* SymPy has functions for calculus: `limit`, `diff`, `integrate`, $\dots$

3
quarto/precalc/logical_expressions.qmd
View File

@@ -269,7 +269,7 @@ A,B = true, false ## also true, true; false, true; and false, false
## Precedence
The question of when parentheses are needed and when they are not is answered by the [precedence](http://julia.readthedocs.org/en/latest/manual/mathematical-operations/#operator-precedence) rules implemented. Earlier, we wrote
The question of when parentheses are needed and when they are not is answered by the [precedence](https://docs.julialang.org/en/v1/manual/mathematical-operations/#Operator-Precedence-and-Associativity) rules implemented. Earlier, we wrote
```{julia}
@@ -625,4 +625,3 @@ In the manual we can read that "In the expression `a && b`, the subexpression `b
answ = 1
radioq(choices, answ)
```

3
quarto/precalc/ranges.qmd
View File

@@ -233,7 +233,7 @@ x = [0, 2, 4, 6, 8, 10]
That is of course only part of the set of even numbers we want. Creating more might be tedious were we to type them all out, as above. In such cases, it is best to *generate* the values.
For this simple case, a range can be used, but more generally a [comprehension](http://julia.readthedocs.org/en/latest/manual/arrays/#comprehensions) provides this ability using a construct that closely mirrors a set definition, such as $\{2k: k=0, \dots, 50\}$. The simplest use of a comprehension takes this form (as we described in the section on vectors):
For this simple case, a range can be used, but more generally a [comprehension](https://docs.julialang.org/en/v1/manual/arrays/#man-comprehensions) provides this ability using a construct that closely mirrors a set definition, such as $\{2k: k=0, \dots, 50\}$. The simplest use of a comprehension takes this form (as we described in the section on vectors):
`[expr for variable in collection]`
@@ -673,4 +673,3 @@ The [product](http://en.wikipedia.org/wiki/Arithmetic_progression) of the terms
val = prod(1:2:19)
numericq(val)
```

5
quarto/precalc/transformations.qmd
View File

@@ -333,13 +333,13 @@ This is off by a fair amount - almost $12$ minutes. Clearly a trigonometric mode
##### Example: a growth model in fisheries
The von Bertalanffy growth [equation](http://www.fao.org/docrep/W5449e/w5449e05.htm) is $L(t) =L_\infty \cdot (1 - e^{k\cdot(t-t_0)})$. This family of functions can be viewed as a transformation of the exponential function $f(t)=e^t$. Part is a scaling and shifting (the $e^{k \cdot (t - t_0)}$) along with some shifting and stretching. The various parameters have physical importance which can be measured: $L_\infty$ is a carrying capacity for the species or organism, and $k$ is a rate of growth. These parameters may be estimated from data by finding the "closest" curve to a given data set.
The von Bertalanffy growth [equation](https://en.wikipedia.org/wiki/Von_Bertalanffy_function) is $L(t) =L_\infty \cdot (1 - e^{k\cdot(t-t_0)})$. This family of functions can be viewed as a transformation of the exponential function $f(t)=e^t$. Part is a scaling and shifting (the $e^{k \cdot (t - t_0)}$) along with some shifting and stretching. The various parameters have physical importance which can be measured: $L_\infty$ is a carrying capacity for the species or organism, and $k$ is a rate of growth. These parameters may be estimated from data by finding the "closest" curve to a given data set.
##### Example: the pipeline operator
In the last example, we described our sequence as scale, over, stretch, and up, but code this in reverse order, as the composition $f \circ g$ is done from right to left. A more convenient notation would be to have syntax that allows the composition of $g$ then $f$ to be written $x \rightarrow g \rightarrow f$. `Julia` provides the [pipeline](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.|>) operator for chaining function calls together.
In the last example, we described our sequence as scale, over, stretch, and up, but code this in reverse order, as the composition $f \circ g$ is done from right to left. A more convenient notation would be to have syntax that allows the composition of $g$ then $f$ to be written $x \rightarrow g \rightarrow f$. `Julia` provides the [pipeline](https://docs.julialang.org/en/v1/manual/functions/#Function-composition-and-piping) operator for chaining function calls together.
For example, if $g(x) = \sqrt{x}$ and $f(x) =\sin(x)$ we could call $f(g(x))$ through:
@@ -645,4 +645,3 @@ q"S(D(f))(n) = f(n) - f(0)"
answ = 1
radioq(choices, answ, keep_order=true)
```

3
quarto/precalc/vectors.qmd
View File

@@ -407,7 +407,7 @@ The last value of a vector is usually denoted by $v_n$. In `Julia`, the `length`
:::{.callout-note}
## More on indexing
There is [much more](http://julia.readthedocs.org/en/latest/manual/arrays/#indexing) to indexing than just indexing by a single integer value. For example, the following can be used for indexing:
There is [much more](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing) to indexing than just indexing by a single integer value. For example, the following can be used for indexing:
* a scalar integer (as seen)
* a range
@@ -990,4 +990,3 @@ q"zs^(1./2)"
answ = 2
radioq(choices, answ, keep_order=true)
```