Nonlinear differential equation struggling to solve

Specific Domains Modelling & Simulations
1

Hi Guys

Please refer to the attached code. The code solves a set of coupled differential equations. The code seems to work well for certain input parameters, however sometimes the system seems to be very stiff and does not solve correctly.

I come from a Matlab background, can anyone provide some guidance on how I can alter setting to get better results?

using DifferentialEquations
using Plots


function statespace!(dx,x,p,t)
    l,m,g,I,kc,kb,C1,C2 = p

    A = x[2]^2*l*sin(x[1])+(-kb*x[3]-C1*x[4])/m
    B = -kc*x[1]-m*g*l*sin(x[1])-C2*x[2]
    C = (1-(m*l^2*cos(x[1])^2)/(l^2*m+I))

    theta_dotdot = ((-m*A*l*cos(x[1]+B)/(l^2*m+I)))*1/C

    x_dotdot = x[2]^2*l*sin(x[1])-theta_dotdot*l*cos(x[1])+(-kb*x[3]-C1*x[4])/m

    dx[1] = x[2]
    dx[2] = theta_dotdot
    dx[3] = x[4]
    dx[4] = x_dotdot  

end

m = 230
l = 5
g = 9.81
I = 230*(l/2)^2
kc = 795774.71
kb = 236*1000
D1 = 0.04*(m+kb)*0
D2 = 0.04*(m+kc)*0

x0 = [0 , 0.0, 0.0, 2.5]  
tspan = (0.0,1)  # time span

p = (m,l,g,I,kc,kb,D1,D2)

prob = ODEProblem(statespace!, x0, tspan, p)
sol = solve(prob,Rodas5(),saveat = 0.001,reltol=1e-8,abstol=1e-8)

plt_1 = plot()
plot!(plt_1, sol.t,sol[1,:],xlabel = "Time (sec)",ylabel = "Angular displacement")
display(plot(plt_1))

plt_2 = plot()
plot!(plt_2, sol.t,sol[2,:],xlabel = "Time (sec)",ylabel = "Angular velocity")
display(plot(plt_2))

plt_3 = plot()
plot!(plt_3, sol.t,sol[3,:],xlabel = "Time (sec)",ylabel = "Linear displacement")
display(plot(plt_3))

plt_4 = plot()
plot!(plt_4, sol.t,sol[4,:],xlabel = "Time (sec)",ylabel = "Linear velocity")
display(plot(plt_4))
2

The simplest thing to try is to experiment with solvers, e.g., a high-order one for stiff systems like RadauIIA9():

sol = solve(prob, RadauIIA9(), reltol=1e-8, abstol=1e-8)

You may also try FBDF, QNDF and the other Rodas solvers like Rodas5P()

2 Likes