CalculusWithJuliaNotes.jl/quarto/ODEs/solve.qmd

# The problem-algorithm-solve interface


{{< include ../_common_code.qmd >}}

This section uses these add-on packages:


```{julia}
using Plots
using MonteCarloMeasurements
```


---


The [DifferentialEquations.jl](https://github.com/SciML) package is an entry point to a suite of `Julia` packages for numerically solving differential equations in `Julia` and other languages. A common interface is implemented that flexibly adjusts to the many different problems and algorithms covered by this suite of packages.

In this section, we review a very informative [post](https://discourse.julialang.org/t/function-depending-on-the-global-variable-inside-module/64322/10) by discourse user `@genkuroki` which very nicely demonstrates the usefulness of the problem-algorithm-solve approach used with `DifferentialEquations.jl`. We slightly modify the presentation below for our needs, but suggest a perusal of the original post.


##### Example: Free fall


The motion of an object under a uniform gravitational field is of interest.


The parameters that govern the equation of motions are the gravitational constant, `g`; the initial height, `y0`; and the initial velocity, `v0`. The time span for which a solution is sought is `tspan`.


A problem consists of these parameters. Typical `Julia` usage would be to create a structure to hold the parameters, which may be done as follows:


```{julia}
struct Problem{G, Y0, V0, TS}
    g::G
    y0::Y0
    v0::V0
    tspan::TS
end

Problem(;g=9.80665, y0=0.0, v0=30.0, tspan=(0.0,8.0)) = Problem(g, y0, v0, tspan)
```

The above creates a type, `Problem`, *and* a default constructor with default values. (The original uses a more sophisticated setup that allows the two things above to be combined.)

Types, as used above, serve two purposes: they bundle together the parameters for later reference and they can be used to dispatch varying methods to solve problems. The `solve` generic in the `Julia` ecosystem dispatches on the type of problem it is given.

In the above code, just calling `Problem()` will create a problem suitable for the earth, passing different values for `g` would be possible for other planets.


To solve differential equations there are many different possible algorithms. Here is the construction of two types to indicate two algorithms:


```{julia}
struct EulerMethod{T}
  dt::T
end
EulerMethod(; dt=0.1) = EulerMethod(dt)

struct ExactFormula{T}
  dt::T
end
ExactFormula(; dt=0.1) = ExactFormula(dt)
```

The above just specifies a type for dispatch-–-the directions indicating what code to use to solve the problem and default constructors. As seen, each constructor specifies a default size for a time step of `0.1`.


A type for solutions is useful for different `show` methods or other methods. One can be created through:


```{julia}
struct Solution{Y, V, T, P<:Problem, A}
    y::Y
    v::V
    t::T
    prob::P
    alg::A
end
```

The different algorithms then can be implemented as part of a generic `solve` function. Following the post we have:


```{julia}
solve(prob::Problem) = solve(prob, default_algorithm(prob))
default_algorithm(prob::Problem) = EulerMethod()

function solve(prob::Problem, alg::ExactFormula)

    (; g, y0, v0, tspan) = prob # property destructuring
    dt = alg.dt                 # direct property access
    t0, t1 = tspan

    ts = range(t0, t1 + dt/2; step = dt)

    y(t) = y0 + v0*(t - t0) - g*(t - t0)^2/2
    v(t) = v0 - g*(t - t0)

    Solution(y.(ts), v.(ts), ts, prob, alg)
end
```

The exact formulas:

$$
\begin{align*}
y(t) &= y_0 + v_0\cdot(t - t_0) - g\cdot(t - t_0)^2/2\\
v(t) &= v_0 - g\cdot(t - t_0),
\end{align*}
$$

are well-known physics formulas, discussed previously, for motion under a constant acceleration. The `ExactFormula` code broadcasts these functions over a range of values in `ts` and then wraps the output up in a `Solution` object so that the answers found can be easily extracted in a uniform manner.

For the Euler method, a `for` loop is utilized to step through the algorithm, in preparation, the new command `fill(y0, n)` is technical. It sets up a storage vector of length `n` which is initially filled with `y0` but for which the second through last are overwritten. There are many other means to do a similar task, including creating an uninitialized vector with `Vector{typeof(y0)}(undef, n)` for which all entries would be subsequently filled in.

```{julia}
function solve(prob::Problem, alg::EulerMethod)

    (; g, y0, v0, tspan) = prob
    dt = alg.dt
    t0, t1 = tspan

    ts = range(t0, t1 + dt/2; step = dt)
    n = length(ts)

    ys = fill(y0, n)
    vs = fill(v0, n)

    for i in 1:n-1
        vs[i+1] = vs[i] - g*dt      # F*h step of Euler
        ys[i+1] = ys[i] + vs[i]*dt  # F*h step of Euler
    end

    Solution(ys, vs, ts, prob, alg)
end
```

Plots of solutions generated by the default values for each method are produced in @fig-projectile-motion-on-the-earth-dt-default.

::: {#fig-projectile-motion-on-the-earth-dt-default}
```{julia}
earth = Problem()
sol_euler = solve(earth)
sol_exact = solve(earth, ExactFormula())

plot(sol_euler.t, sol_euler.y;
     label="Euler's method (dt = $(sol_euler.alg.dt))", linestyle=:auto)
plot!(sol_exact.t, sol_exact.y;
      label="exact solution", linestyle=:auto)

title!("On the Earth"; xlabel="t", legend=:bottomleft)
```

Plot of exact and approximate solutions, the latter using the default time step size
:::

Following the post, since the time step `dt = 0.1` is not small enough, the error of the Euler method is readily identified in @fig-projectile-motion-on-the-earth-dt-default.

Next we change the algorithm's default parameter for `dt` to be smaller. @fig-projectile-motion-on-the-earth-dt-modified shows a much improved agreement between the exact answer and the approximate one found with `EulerMethod`.

::: {#fig-projectile-motion-on-the-earth-dt-modified}
```{julia}
earth₂ = Problem()
sol_euler₂ = solve(earth₂, EulerMethod(dt = 0.01))
sol_exact₂ = solve(earth₂, ExactFormula())

plot(sol_euler₂.t, sol_euler₂.y;
     label="Euler's method (dt = $(sol_euler₂.alg.dt))", linestyle=:auto)
plot!(sol_exact₂.t, sol_exact₂.y;
      label="exact solution", linestyle=:auto)

title!("On the Earth"; xlabel="t", legend=:bottomleft)
```

Model of projectile motion on the earth along with exact solution. This approximation used a modification from the default for the step size.
:::

The code is mostly a template. It is worth noting that only one line of code was modified, and in that line only the method required a modification.


Were the moon to be considered, the gravitational constant would need adjustment. This parameter is a property of the problem, not the solution algorithm, as `dt` is.


Such adjustments are made by passing different values to the `Problem` constructor. Again, just the one line needs modification.

::: {#fig-projectile-motion-on-the-moon}
```{julia}
moon = Problem(g = 1.62, tspan = (0.0, 40.0))

sol_eulerₘ = solve(moon)
sol_exactₘ = solve(moon, ExactFormula(dt = sol_euler.alg.dt))

plot(sol_eulerₘ.t, sol_eulerₘ.y;
     label="Euler's method (dt = $(sol_eulerₘ.alg.dt))", linestyle=:auto)
plot!(sol_exactₘ.t, sol_exactₘ.y;
      label="exact solution", linestyle=:auto)

title!("On the Moon"; xlabel="t", legend=:bottomleft)
```

Model for projectile motion on the moon. Only modest changes needed to be introduced from the model for the motion on oearth illustrated in @fig-projectile-motion-on-the-earth-dt-default.
:::

The code above also adjusts the time span in addition to the graviational constant. The algorithm for the exact formula is set to use the `dt` value used in the `euler` formula, for easier comparison. Otherwise, outside of the labels, the patterns are the same. Only those things that need changing are changed, the rest comes from defaults.


The above shows the benefits of using a common interface---new problems can be approached through only minor adjustments to the parameters, yet the calling pattern remains the same.

Next, the post illustrates how *other* authors *could* extend this code. The `solve` method dispatches on the problem type and the method type. Adding a new method to `solve` requires defining new method type and the algorithm for that type in the extension of `solve`.

For example, the following adds a sympletic method which conserves a quantity, allowing the approximate solutions to track long-term evolution without drift.


```{julia}
struct Symplectic2ndOrder{T}
    dt::T
end
Symplectic2ndOrder(; dt=0.1) = Symplectic2ndOrder(dt)

function solve(prob::Problem, alg::Symplectic2ndOrder)

    g, y0, v0, tspan = prob.g, prob.y0, prob.v0, prob.tspan
    dt = alg.dt
    t0, t1 = tspan

    ts = range(t0, t1 + dt/2; step = dt)
    n = length(ts)

    ys = fill(y0, n)
    vs = fill(v0, n)

    for i in 1:n-1
        vs[i+1] = vs[i] - g*dt
        ys[i+1] = ys[i] + (vs[i] + vs[i+1])/2 * dt
    end

    Solution(ys, vs, ts, prob, alg)
end
```

Had the two prior methods been in a package, the other user could still extend the interface, as above, with just a slight standard modification.

The exact same approach to solving a problem works for this new type:

::: {#fig-projectile-motion-symplectic-2-order}
```{julia}
earth₃ = Problem()
sol_sympl₃ = solve(earth₃, Symplectic2ndOrder(dt = 2.0))
sol_exact₃ = solve(earth₃, ExactFormula())

plot(sol_sympl₃.t, sol_sympl₃.y;
     label="2nd order symplectic (dt = $(sol_sympl₃.alg.dt))", linestyle=:auto)
plot!(sol_exact₃.t, sol_exact₃.y;
      label="exact solution", linestyle=:auto)

title!("On the Earth"; xlabel="t", legend=:bottomleft)
```

Plot of exact solution and approximate solution using a $2$nd-order symplectic method
:::

Finally, the author of the post shows how the interface can compose with other packages in the `Julia` package ecosystem. This example uses the external package `MonteCarloMeasurements` which plots the behavior of the system for perturbations of the initial value, as seen in @fig-illustration-with-monte-carlo-measurements.

::: {#fig-illustration-with-monte-carlo-measurements}
```{julia}
using MonteCarloMeasurements # introduces ± operation

earth₄ = Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)
sol_euler₄ = solve(earth₄)
sol_sympl₄ = solve(earth₄, Symplectic2ndOrder(dt = 2.0))
sol_exact₄ = solve(earth₄, ExactFormula())


P = plot(sol_euler₄.t, sol_euler₄.y;
         label="Euler's method (dt = $(sol_euler₄.alg.dt))", linestyle=:auto)

Q = plot(sol_sympl₄.t, sol_sympl₄.y;
         label="2nd order symplectic (dt = $(sol_sympl₄.alg.dt))", linestyle=:auto)

R = plot(sol_exact₄.t, sol_exact₄.y;
         label="exact solution", linestyle=:auto)

title!.((P,Q,R), "On the Earth"; xlabel="t", legend=:bottomleft, ylims=(-100, 60))

plot(P, Q, R)
```

Figures showing solutions to a differential equation where the initial values have some specified uncertainty
:::

The only change was in the problem, `Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)`, where a different number type is used which accounts for uncertainty. The rest follows the same pattern.


This example, shows the flexibility of the problem-algorithm-solve pattern while maintaining a consistent pattern for execution.