# 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

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)

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

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.

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

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.

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)

return sol
end

sol = run_example()

times = LinRange(0.0, 10.0, 100)

z = [sol(t)[3] for t in times]
w = [sol(t)[6] 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:

(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!

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.