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?
Consider:
\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)
where
\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$")
plot!(xlabel=L"t")