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).
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 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.
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).
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 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.
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.
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:
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.
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)
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:
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:
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:
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:
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]`:
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`.
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.)
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 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.
(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:
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.
approximations as well as the gradient." On this problem it performs similarly to `Newton`, though in general may be preferable for higher-dimensional problems.
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 `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:
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.
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]`).
\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``:
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``.
