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

WIP

Browse Source
This commit is contained in:
jverzani
2022-08-21 10:15:34 -04:00
parent 65db6baa3a
commit e901163fcf
8 changed files with 645 additions and 496 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -0,0 +1,292 @@
# The SciML suite of packages
The `Julia` ecosystem advances rapidly and 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. Many recent efforts of this organization have been to write uniform interfaces to other packages in the ecosystem.
The basic interface of many of 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).
## Symbolic math (`Symbolics`)
The `Symbolics` package, along with `SymbolicUtils` and `ModelingToolkit` are provided by this organization. See this section on [Symbolics](./symbolics.html) for additional details or the package [documentation](https://symbolics.juliasymbolics.org/stable/) or the documentation for [SymbolicsUtils](https://github.com/JuliaSymbolics/SymbolicUtils.jl).
## Solving equations (`LinearSolve`, `NonlinearSolve`)
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.
The `NonlinearSolve` package can be seen as an alternative to our use of the `Roots` package. 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:
```julia
using NonlinearSolve
```
Unlike `Roots` the package handles problems beyond the univariate case, as such the simplest problems have a little extra 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 use the "dots" to allow vectorization of the basic math operations, as `u` is a container of values. The `@.` macro makes this quite easy, as illustrated above.
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 done automatically.
```julia
soln = solve(prob, NewtonRaphson())
```
The basic interface for retrieving the solution from the solution object is to use indexing:
```julia
soln[]
```
----
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())
```
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``. 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...)
```
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:
```julia
λ = (x, p=nothing) -> NonlinearSolve.ForwardDiff.gradient(peaks, x)
u0 = @SVector[1.0, 1.0]
prob = NonlinearProblem(λ, u0)
u = solve(prob, NewtonRaphson())
```
We can see that this value is a "zero" through:
```julia; error=true
λ(u.u)
```
### Using Modeling toolkit to model the non-linear problem
Nonlinear problems can also be approached symbolically, using `ModelingToolkit`. 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. Here 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. But 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 straight. The `@named` macro handles this behind the scenes by simply rewriting the syntax:
```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 `a => 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. The two packages mentioned here are `Optim` and `BlackBoxOptimization`. The tie-in packages `OptimizationOptimJL` and `OptimizationBBO` must be loaded for these.
### Local
We begin with a simple example from first semester calculus: among all rectangle of perimeter 25, find the one with the largest area. The mathematical setup has a constraint (``P=25=2x+2y``) and from the objective (``A=xy``), 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
A(x, p) = @.(- x * (25 - 2x)/2)
```
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 aim to use:
```julia
using Optimization
using OptimizationOptimJL
```
Next, we define an optimization function with information on how its derivatives will be take. The following uses `ForwardDiff`, which is a good choice in the typical calculus setting with a small number of inputs (just ``q`` 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 form for 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 `minimum` property holds the identified minimum
```julia
soln.minimum
```
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 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, which are needed in the choice of method, `Newton`, that follows.
Solving this problem the 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, solving the the minimum perimeter rectangle for a fixed area (``25`` below), 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, BFGS())
```
We used an initial guess of ``x=4`` above. The `BFGS` method is
described in the documentation as "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`.
### Two dimensional
## Integration (`Integrals.jl`)