Apologies if this question is poorly worded, I’m very new to numerical continuation and nonlinear frequency analysis. I have a nonlinear system of ODEs, and I’d like to numerically find periodic solutions for a changing parameter. Ideally, I’d like to recreate Figure 1 in \verb|[1]| for the Duffing equation. I haven’t included that figure here, I don’t think it’s necessary. Basically, given a nonlinear system, I want to find a periodic (steady state) response for a range of (single) parameter inputs. Looking at the Duffing Equation, which is basically a forced spring-mass damper with a nonlinearity thrown in, I’d like to keep all parameters constant except for \omega, and then numerically find a one-period-long time-domain solution to the nonlinear dynamics for each value of \omega.
Duffing Equation
\ddot{z} + c \dot{z} + k z + \alpha z^3 = A \cos \left(\omega t \right)
Duffing Equation (redefined)
- We can get rid of t by rewriting the dynamics as a fourth-order autonomous system
\begin{align} \dot{z}_1 &= z_2 \\ \dot{z}_2 &= -k z_1 - c z_2 -\alpha z_1^3 + A z_4 \\ \dot{z}_3 &= z_3 + \omega z_4 - z_3 \left(z_3^2 + z_4^2 \right) \\ \dot{z}_4 &= -\omega z_3 + z_4 - z_4 \left( z_3^2 + z_4^2 \right) \end{align}
Note that the solution of z_4 is equivalent to \cos (\omega t), as shown in the \dot{z_2} expression.
For many values of \omega (and all other parameters held constant), I would like to find a periodic solution for the steady state response of the system. I’m trying to accomplish this with the code below.
using Setfield
using LinearAlgebra
using BifurcationKit
using ModelingToolkit
using DifferentialEquations
@variables t z[1:4](t)
@parameters c k α A ω
@derivatives D'~t
duffing, f, J = let
eqs = [
D(z[1]) ~ z[2],
D(z[2]) ~ -k*z[1] - c*z[2] - α*z[1]^3 + A*z[4],
D(z[3]) ~ z[3] + ω*z[4] - z[3]*(z[3]^2 + z[4]^2),
D(z[4]) ~ -ω*z[3] + z[4] - z[4] * (z[3]^2 + z[4]^2)
]
# Build ModelingToolkit system
duffing = ODESystem(eqs, t, z, [ω, A, c, k, α]; name=:Duffing)
# Build function for dynamics
func = eval(first(generate_function(duffing)))
f = (x,u,t=0)->func(x, (u), t)
# Build function for Jacobian
jac = eval(first(generate_jacobian(duffing)))
J = (x,u,t=0)->jac(x,(u), t)
duffing, f, J
end
# Time slices
t = collect(0.0:0.01:10.0)
# Guess for orbit period
T = 1.0
# Guess for periodic orbit
orbitguess = vcat(sin.(t), cos.(t), ones(length(t) * 2), T)
# Construct Periodic Orbit Trapezoidal Problem
pb = PeriodicOrbitTrapProblem(f, J, 0.0, 0.0, 10)
newton(pb, orbitguess, (2.5, 2.5, 0.2, 1.0, 0.05), NewtonPar())
This code errors on the newton call, because pb seems to think the number of states, N, is zero. Am I setting up PeriodicOrbitTrapProblem incorrectly?
Thanks in advance for the help.
And thanks for the package @rveltz!
Reference
\verb|[1]| Nguyen, Duc H., Mark H. Lowenberg, and Simon Neild. “Nonlinear Frequency Response Analysis to Inform Aircraft Control Law Design.” In AIAA Scitech 2020 Forum, p. 0605. 2020.
(but you had no other choice for now)
