Newmark beta method

[termi-official]

Prototype to resolve #411 .

Checklist

• [ ] Appropriate tests were added
• [ ] Any code changes were done in a way that does not break public API
• [ ] All documentation related to code changes were updated
• [ ] The new code follows the contributor guidelines, in particular the SciML Style Guide and COLPRAC.
• [ ] Any new documentation only uses public API

[termi-official]

The last remaining problem is a bit tricky. For a given DynamicalODEFunction Newmark needs to solve for a "future accelleration" by solving the problem $Δtₙ f₁(innertmp₂ + z β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z$ for $z$ with given $innertmp, \alpha,\beta$.

The first issue I run into is passing down the analytical Jacobian:

using OrdinaryDiffEq
function f1_harmonic(dv, v, u, p, t)
    dv .= -u
end
function f2_harmonic(du, v, u, p, t)
    du .= v
end
harmonic_jac_f1(J, v, u, p, t) = J[1,1] = -1.0
ff_harmonic = DynamicalODEFunction(f1_harmonic, f2_harmonic; jac = harmonic_jac_f1)
DiffEqBase.has_jac(ff_harmonic) # false
DiffEqBase.has_jac(ff_harmonic.f1) # false
ff_f1 = ODEFunction{false}(ODEFunction{false}(f1_harmonic); jac=harmonic_jac_f1!)
ff_harmonic = DynamicalODEFunction(ff_f1, f2_harmonic; jac = harmonic_jac_f1)
DiffEqBase.has_jac(ff_harmonic) # false
DiffEqBase.has_jac(ff_harmonic.f1) # false 

The second problem is getting the AD on track. My current strategy is to pass a helper type ArrayPartitionNLSolveHelper around, which constructs an ArrayPartition when multiplied with a Vector. The idea is to provide an interface for inplace evaluation of f.f1 (DynamicalODEFunction) which is compatible with the evaluation form required by the internal nonlinear solver, such that we can call $f₁(innertmp + z \alpha,t)$, where $\alpha$ is the mentioned ArrayPartitionNLSolveHelper. Should we head down this road or is there possibly a better design?

[termi-official]

I think I can call this now a first prototype. I am not very sure why it manages to converge to the correct solution, but I assume the residual is just wrong up to scaling. Another problem is that the Jacobian is per construction wrong. The derivative of $f_1$ w.r.t. the input parameter is a rectangular matrix. For Newmark we assume that the input space is partitioned into 2 equally sized parts (velocity $v$ and displacement $u$) and that $f_2$ is just $du = v$. The first assumption tells us that the Jacobian is blocked with 2 blocks. For the inner nonlinear solve in Newmark we now obtain the system matrix for the Newton as $\partial_z f_1(u+\gamma_1 z,v+\gamma_2 z) = \gamma_1 J_{1,1}(u+\gamma_1 z,v+\gamma_2 z) + \gamma_2 J_{1,2}(u+\gamma_1 z,v+\gamma_2 z)$ where $J_{i,j}$ is the i,j-th block of the Jacobian of $f$.

However, this seems to break a few assumptions taken in the code. Hence I could need some feedback on how to proceed before putting time into a design which we do not want in the library. I already think that ArrayPartitionNLSolveHelper is not great, but I cannot figure out a better design by myself either.

Furthermore I cannot figure out how to integrate the W-transformation here. Any pointers?