ODE solution in Symbolics/SymbolicNumericIntegration vs. ModelingToolkit

I’m comparing the analytic and numeric solutions of a linear, second order ODE, and get some weird (arbitrary) result with the analytic solution found by using SymbolicNumericIntegration… I assume I’m using the tools incorrectly…

Here is the ODE I consider:

With x_1(0) = 1 and x_2(0) = -1, and a unit step input for u(t) at t=0, the solutions I find are:

where superscript “a” indicates analytic solution. NOTE that the analytic (solid) and numeric (dashed) solutions seem to diverge at the end of the time period.

The system is unstable, so a divergence might be possible. However, every time I plot the analytic solution, I get a different result. I’m not talking about re-solving the solution analytically; I’m talking about re-plotting the found solution. So this indicates that there is something wrong with the way I represent the analytic solution.

  • What is it I do that is wrong?


\frac{dx}{dt} = Ax+Bu

leading to:

x(t) = \Phi(t)\cdot \left( x(0) + \int_0^t \Phi(-\tau)Bu(\tau)d\tau \right)


\Phi(t) = \exp(At)

My code is as follows:

using SymbolicNumericIntegration
using Plots

@variables τ t

Φ(t) = [1//4*exp(-10t) + 3//4*exp(t) -3//4*exp(-10t) + 3//4*exp(t);
        -1//4*exp(-10t) + 1//4*exp(t) 3//4*exp(-10t) + 1//4*exp(t)]
B = [1;0]
x0 = [1; -1]

# constant input
intg(τ) = [i[1] for i in integrate.(Φ(-τ)*B,τ)]  # stripping off errors, etc.

intg_def(t) = intg(t) .- substitute.(intg(τ),(Dict(τ=>0),))

x(t) = simplify.(expand.(Φ(t)*x0 .+ Φ(t)*intg_def(t)))

Tf = 3.8
plot(x(t)[1], xlim=(0,Tf),lw=2.5, lc=:blue,la=0.3,label=L"x_1^\mathrm{a}")
plot!(x(t)[2], xlim=(0,Tf),lw=2.5, lc=:red,la=0.3,label=L"x_2^\mathrm{a}")
plot!(title=L"Diagonalizable system with input $u(t)=1$")

Always check the analytical result by using differentiation.

(We should do that automatically and throw a warning)

It is not completely clear to me how to do the differentiation… my solution x(t) is a function of symbolic variable t. Doing Differentiate(x) of x(t) produces 0 (zero), while Differentiate(x(t)) gives an error message…

Because the analytic result is almost identical to the numeric result over most of the time span, and because if I re-do the plotting of the analytic results, there is only a slight “wiggling” at the largest values of the time span, I think the solution is correct… so I’m suspecting something wrong in the way I have created the functions – or possibly that the plotting recipe doesn’t understand the function properly.