# Converting a toy example to DifferentialEquations

**URL:** <https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [May 13, 2021, 8:44am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070 "2021-05-13T08:44:19Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 8:44am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/1 "2021-05-13T08:44:20Z")

</div>

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

```nohighlight
"""
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](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.

---

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 13, 2021, 9:10am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/2 "2021-05-13T09:10:31Z")

</div>

I think [this example](https://diffeq.sciml.ai/stable/tutorials/dae_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](https://diffeq.sciml.ai/stable/solvers/dae_solve/#ODEInterfaceDiffEq.jl)

---

<div class="post-metadata">

**Author:** ![maajdl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maajdl/32/22838_2.png) [@maajdl](https://discourse.julialang.org/u/maajdl)\
**Post date:** [May 13, 2021, 9:14am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/3 "2021-05-13T09:14:35Z")

</div>

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

[https://diffeq.sciml.ai/stable/solvers/dae\_solve/](https://diffeq.sciml.ai/stable/solvers/dae_solve/)

> #### [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.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 9:16am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/4 "2021-05-13T09:16:47Z")

</div>

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: [[1406.6218] Dynamic Model of a Pumping Kite Power System](http://arxiv.org/abs/1406.6218)

---

<div class="post-metadata">

**Author:** ![maajdl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maajdl/32/22838_2.png) [@maajdl](https://discourse.julialang.org/u/maajdl)\
**Post date:** [May 13, 2021, 9:19am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/5 "2021-05-13T09:19:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 9:29am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/6 "2021-05-13T09:29:08Z")

</div>

This is the Julia code I have so far: [KiteViewer/KPS3.jl at sim · ufechner7/KiteViewer · GitHub](https://github.com/ufechner7/KiteViewer/blob/sim/src/KPS3.jl)

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

```julia
        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… 🙂

---

<div class="post-metadata">

**Author:** ![maajdl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maajdl/32/22838_2.png) [@maajdl](https://discourse.julialang.org/u/maajdl)\
**Post date:** [May 13, 2021, 9:45am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/7 "2021-05-13T09:45:04Z")

</div>

Could this help?

[https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/src/dae\_premade\_problems.jl](https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/src/dae_premade_problems.jl)

Found that here:

[https://diffeq.sciml.ai/stable/types/dae\_types/#DiffEqProblemLibrary.DAEProblemLibrary.prob\_dae\_resrob](https://diffeq.sciml.ai/stable/types/dae_types/#DiffEqProblemLibrary.DAEProblemLibrary.prob_dae_resrob)

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 13, 2021, 9:59am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/8 "2021-05-13T09:59:07Z")

</div>

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 (??).

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 10:13am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/9 "2021-05-13T10:13:13Z")

</div>

Thank you, I will have a look.

I converted my toy example:

```julia
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, :]

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 10:30am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/10 "2021-05-13T10:30:13Z")

</div>

How could I rewrite this example using a mass matrix?

---

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 13, 2021, 10:59am UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/11 "2021-05-13T10:59:00Z")

</div>

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

```julia
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)[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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 12:13pm UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/12 "2021-05-13T12:13:34Z")

</div>

> [@ufechner7](#):
>
> 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?

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](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/](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/](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.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 12:50pm UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/13 "2021-05-13T12:50:49Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 13, 2021, 1:14pm UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/14 "2021-05-13T13:14:52Z")

</div>

You can control that with `saveat`.  
[https://diffeq.sciml.ai/stable/basics/common\_solver\_opts/#Output-Control](https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Output-Control)

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 13, 2021, 1:40pm UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/15 "2021-05-13T13:40:59Z")

</div>

Now also plotting works nicely:

```julia
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()

```

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [May 13, 2021, 3:45pm UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/16 "2021-05-13T15:45:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 5:16pm UTC](https://discourse.julialang.org/t/converting-a-toy-example-to-differentialequations/61070/17 "2021-05-13T17:16:23Z")

</div>

> [@tomaklutfu](#):
>
> 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.