 # Converting a toy example to DifferentialEquations

I am trying to convert the following Python example to Julia:

"""
Tutorial example showing how to use the implicit solver RADAU.
It simulates a falling mass.
"""
import numpy as np
import pylab as plt
from assimulo.problem import Implicit_Problem # Imports the problem formulation from Assimulo
from assimulo.solvers import Radau5DAE        # Imports the solver RADAU from Assimulo

G_EARTH  = np.array([0.0, 0.0, -9.81]) # gravitational acceleration

# Example one:
# Falling mass.
# State vector y   = mass.pos, mass.vel
# Derivative   yd  = mass.vel, mass.acc
# Residual     res = (y.vel - yd.vel), (yd.acc - G_EARTH)
def res1(t, y, yd):
res_0 = y[3:6]  - yd[0:3]
res_1 = yd[3:6] - G_EARTH
return np.append(res_0, res_1)

def run_example():
# Set the initial conditons
t0  = 0.0                   # Initial time
vel_0 = [0.0, 0.0, 50.0]    # Initial velocity
pos_0 = [0.0, 0.0,  0.0]    # Initial position
acc_0 = [0.0, 0.0, -9.81]   # Initial acceleration
y0  = pos_0 + vel_0         # Initial pos, vel
yd0 = vel_0 + acc_0         # Initial vel, acc

model = Implicit_Problem(res1, y0, yd0, t0) # Create an Assimulo problem
model.name = 'Falling mass' # Specifies the name of problem (optional)

sim = Radau5DAE(model)      # Create the Radau solver

tfinal = 10.0           # Specify the final time
ncp    = 500            # Number of communcation points (number of return points)

# Use the .simulate method to simulate and provide the final time and ncp (optional)
time, y, yd = sim.simulate(tfinal, ncp)

# plot the result
pos_z = y[:,2]
vel_z = y[:,5]


The Python package Assimulo has a very good documentation for solving implicit problems:
https://jmodelica.org/assimulo/tutorial_imp.html

I did not found anything similar in the documentation of DifferentialEquations.

Which problem shall I define to solve this example?
Shall I create a DAEProblem? But then I cannot use Radau, correct?

Please keep in mind that my real problem is large, I coded it in Numba/Python (now converting it to Julia) and would prefer to stick to the Radau solver, because it performs well for my stiff problem of about 50 equations.

I think this example is analogous to your link.

You can specify the problem as an ODEProblem with a mass matrix if you want to use Radau.
https://diffeq.sciml.ai/stable/solvers/dae_solve/#ODEInterfaceDiffEq.jl

Have you checked this page and the entry “FIRK-Methods”?

https://diffeq.sciml.ai/stable/solvers/dae_solve/

#### FIRK Methods

• RadauIIA5 - An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.

To be honest, I have no idea what a “mass matrix” is. Is it always possible to convert an implicit problem into mass matrix form?

You can find the problem that I am trying to implement here: http://arxiv.org/abs/1406.6218

The Python code might be easier to read than the paper ! It’s not immediate to see where the “implicit” comes from.
And why a DAE solver is needed or if a “simpler” ODE solver would be enough.

This is the Julia code I have so far: KiteViewer/KPS3.jl at sim · ufechner7/KiteViewer · GitHub

It calculates the residual based on the vectors y and yd.

        N-point tether model:
Inputs:
State vector state_y   = pos0, pos1, ..., posn-1, vel0, vel1, ..., veln-1
Derivative   der_yd    = vel0, vel1, ..., veln-1, acc0, acc1, ..., accn-1
Output:
Residual     res = res0, res1 = pos0,  ..., vel0, ...


And I now have to find a way to connect this code (actually the function loop()) to a solver.

But I want to start with the toy example… In MATLAB (and Julia?), a “mass matrix” is a matrix M in an implicit differential equation

M(x;t)\frac{\mathrm{d}x}{\mathrm{d}t} = f(x,u;t)

When M is rank deficient, this is effectively a DAE. Originally, this was the way to specify DAEs in MATLAB.

Sometimes, such models can be written in semi-explicit form as

\frac{\mathrm{d}x_1}{\mathrm{d}t} = f_1(x_1,x_2,u;t) \\ 0 = f_2(x_1,x_2,u;t)

which may be simpler to write down. But perhaps the version with mass matrix is more efficient to solve (??).

Thank you, I will have a look.

I converted my toy example:

using DifferentialEquations, Sundials
# Tutorial example showing how to use an implicit solver
# It simulates a falling mass.

const G_EARTH  = [0.0, 0.0, -9.81]

# Falling mass.
# State vector y   = mass.pos, mass.vel
# Derivative   yd  = mass.vel, mass.acc
# Residual     res = (y.vel - yd.vel), (yd.acc - G_EARTH)
function res1(res, yd, y, p, t)
res[1:3] .= y[4:6] - yd[1:3]
res[4:6] .= yd[4:6] - G_EARTH
end

vel_0 = [0.0, 0.0, 50.0]    # Initial velocity
pos_0 = [0.0, 0.0,  0.0]    # Initial position
acc_0 = [0.0, 0.0, -9.81]   # Initial acceleration
y0  = vcat(pos_0, vel_0)    # Initial pos, vel
yd0 = vcat(vel_0, acc_0)    # Initial vel, acc

tspan = (0.0, 10.0)         # time span

differential_vars = [true, true, true, true, true, true]
prob = DAEProblem(res1, yd0, y0, tspan, differential_vars=differential_vars)

y = solve(prob, IDA())

pos_z = y[3, :]
vel_z = y[6, :]


How could I rewrite this example using a mass matrix?

Your toy problem isn’t actually a DAE so it may not be the best example, but this works:

using OrdinaryDiffEq, LinearAlgebra, Plots

function res1!(du, u, p, t)
# p is gravity vector
for i in 1:3
du[i] = u[i+3]
du[i+3] = p[i]
end
end

function run_example()
G_EARTH = [0.0 0.0 -9.81]

tspan = (0.0, 10.0)

vel_0 = [0.0, 0.0, 50.0]    # Initial velocity
pos_0 = [0.0, 0.0,  0.0]    # Initial position
y0    = vcat(pos_0, vel_0)  # Initial pos, vel

M = Matrix(1I, 6, 6)

ode_f = ODEFunction(res1!, mass_matrix = M)
prob  = ODEProblem(ode_f, y0, tspan, G_EARTH)
sol   = solve(prob, RadauIIA5())

return sol
end

sol = run_example()

times = LinRange(0.0, 10.0, 100)

z = [sol(t) for t in times]
w = [sol(t) for t in times]

display(plot(times, z, label = "z"))
display(plot!(times, w, label = "z dot"))


There’s probably a more elegant way to make those plots.

1 Like

Yes, it’s similar to lowering the order of an ODE. You can always define new variables for things like x^2 to expand the system and end up with a linearly-implicit DAE (i.e. a mass matrix ODE where M is singular). You can see how that represents a DAE in this tutorial:

https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Handling-Mass-Matrices

(which we may want to split out to a new tutorial since it’s pretty deep). That’s as opposed to the fully-implicit DAE form:

https://diffeq.sciml.ai/stable/tutorials/dae_example/

Some solvers (BDF) naturally do the latter, while others (implicit RKs) naturally do the former, hence the two forms. The docs will probably get a heavy makeover soon as ModelingToolkit’s symbolic stuff is really changing this landscape (and automates the conversion between these two and a third DAE form). For a bit on that you may want to check out:

https://mtk.sciml.ai/dev/mtkitize_tutorials/modelingtoolkitize_index_reduction/

I’m still a bit uncertain how to weave all of the symbolics into the DiffEq docs in a nice way.

Thank you!

One additional question:
The time vector that is returned in the solution shows very inequal time steps.

Is there a possibility to get an equal-distant time vector of a pre-defined size as result?

You can control that with saveat.
https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Output-Control

2 Likes

Now also plotting works nicely:

using DifferentialEquations, Sundials, GLMakie
# Tutorial example showing how to use an implicit solver
# It simulates a falling mass.

const G_EARTH  = [0.0, 0.0, -9.81]

# Falling mass.
# State vector y   = mass.pos, mass.vel
# Derivative   yd  = mass.vel, mass.acc
# Residual     res = (y.vel - yd.vel), (yd.acc - G_EARTH)
function res1(res, yd, y, p, t)
res[1:3] .= y[4:6] - yd[1:3]
res[4:6] .= yd[4:6] - G_EARTH
end

vel_0 = [0.0, 0.0, 50.0]    # Initial velocity
pos_0 = [0.0, 0.0,  0.0]    # Initial position
acc_0 = [0.0, 0.0, -9.81]   # Initial acceleration
y0  = vcat(pos_0, vel_0)    # Initial pos, vel
yd0 = vcat(vel_0, acc_0)    # Initial vel, acc

tspan = (0.0, 10.2)         # time span

differential_vars = ones(Bool, 6)
prob = DAEProblem(res1, yd0, y0, tspan, differential_vars=differential_vars)

sol = solve(prob, IDA(), saveat=0.1)

time = sol.t
y = sol.u

pos_z = sol[3, :]
vel_z = sol[6, :]

# plot the result
f = Figure()
ax1 = Axis(f[1, 1], yticklabelcolor = :blue, xlabel="time [s]", ylabel = "pos_z [m]")
ax2 = Axis(f[1, 1], yticklabelcolor = :red, yaxisposition = :right, ylabel = "vel_z [m/s]")
lines!(ax1, time, pos_z, color=:green)
lines!(ax2, time, vel_z, color=:red)
current_figure()


You can also get this values from solution directly by call syntax sol(t0:0.1:tend) if you did not specify saveat. However, its support is not available for all integrator methods or as good as others if I understand correctly.

saveat uses the interpolator. Each method’s interpolation has different properties.