Duffing Equation

The source files for all examples can be found in /examples.

The Duffing equation is a nonlinear second-order ordinary differential equation commonly used to model various physical systems exhibiting nonlinear behavior, such as electrical circuits, mechanical systems, and biological systems. From the mechanical domain point of view, it is an extension of the simple mass-spring-damper system that incorporates a cubic spring. Despite its simple form, the Duffing equation exhibits chaotic behavior. The equation is given by:

\[m\,\frac{d^2 x}{dt^2} + c\,\frac{d x}{dt} + k\,x + k_3\,x^3 = F(t),\]

where $x$ represents the position, $m$ is the mass, $c$ is the damping coefficient, $k$ is the spring stiffness, $k_3$ is the cubic spring stiffness, and $F(t)$ denotes an external force applied to the system.

In standard bond graphs, the compliance element is linear, while for the Duffing equation, we need a nonlinear compliance for the cubic spring. Therefore, we need to introduce a nonlinear compliance element to model the Duffing equation using the Bond Graph Toolkit. The library allows users to define custom elements and utilize the Power connector as interface.

Bond Graph Toolkit

First, we import the Bond Graph Toolkit module to model the system and the DifferentialEquations.jl package to solve the resulting ordinary differential equation (ODE). Additionally, we import the independent variable t from the Bond Graph Toolkit, which will be used to define a custom forcing term.

using BondGraph
using BondGraph: t, D
using DifferentialEquations
using ModelingToolkit
using Plots
using Symbolics.Latexify

Building the model

Custom compliance

We now define a custom element called Spring3 to include the non-linearity of the Duffing equation. The element is defined by the function below.

function Spring3(; name, k = 1.0, x = 0.0)
    # Initialize the Power connector
    @named power = Power()

    # Initialize the displacement state with its initial value `x`
    @variables q(t) = x
    # Set the equation parameters
    ps = @parameters C = 1 / k

    # Define Spring3 element equations with the non-lineariyy of the cubic spring
    eqs = [
        power.e ~ q^3 / C
        D(q) ~ power.f
    ]

    # Build the system with the equations, state and connector.
    compose(ODESystem(eqs, t, [q], ps; name = name), power)
end

Spring3 (generic function with 1 method)

The Spring3 function creates a bond graph subsystem representing the cubic spring. You can specify the following parameters:

name: The name of the Spring3 element (default: "").
k: The cubic spring stiffness constant (default: 1.0).
x: The initial displacement of the spring (default: 0.0).

Setting the elements

Now we initialize and define the cubic spring custom element and the remaining single port elements of the system (mass, spring, and damper) as follows:

The parameters values

α = 1.0
β = 5.0
δ = 0.02
γ = 8.0
ω = 0.5

0.5

Define the elements

@named m = Mass(m = 1.0)
@named s = Spring(k = α)
@named s3 = Spring3(k = β)
@named d = Damper(c = δ)

\[ \begin{align} power_{+}e\left( t \right) =& R power_{+}f\left( t \right) \end{align} \]

The forcing term can be set as an expression or even a Julia function. In this example, we will first intialize a generic forcing term that will be defined before solving the ODE.

@parameters F(t)
@named f = Se(F)

\[ \begin{align} power_{+}e\left( t \right) =& F\left( t \right) \end{align} \]

Then, we build the system by simply passing the one-port elements to the 1-junction function as arguments.

@named mdl = Junction1(m, d, s, s3, [-1, f])

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( m_{+}power, power \right) \\ \mathrm{connect}\left( d_{+}power, power \right) \\ \mathrm{connect}\left( s_{+}power, power \right) \\ \mathrm{connect}\left( s3_{+}power, power \right) \\ \mathrm{connect}\left( power, f_{+}power \right) \\ \frac{\mathrm{d} m_{+}power_{+}f\left( t \right)}{\mathrm{d}t} = \frac{m_{+}power_{+}e\left( t \right)}{m_{+}I} \\ d_{+}power_{+}e\left( t \right) = d_{+}R d_{+}power_{+}f\left( t \right) \\ s_{+}power_{+}e\left( t \right) = \frac{s_{+}q\left( t \right)}{s_{+}C} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} = s_{+}power_{+}f\left( t \right) \\ s3_{+}power_{+}e\left( t \right) = \frac{\left( s3_{+}q\left( t \right) \right)^{3}}{s3_{+}C} \\ \frac{\mathrm{d} s3_{+}q\left( t \right)}{\mathrm{d}t} = s3_{+}power_{+}f\left( t \right) \\ f_{+}power_{+}e\left( t \right) = f_{+}F\left( t \right) \\ \end{array} \right] \end{equation} \]

Analysing the model

We can visualize the bond-graph connections by generating a graph plot of the model.

generate_graph(mdl)

We can obtain the bond-graph equations using the code below.

equations(expand_connections(mdl))

\[ \begin{align} \frac{\mathrm{d} m_{+}power_{+}f\left( t \right)}{\mathrm{d}t} =& \frac{m_{+}power_{+}e\left( t \right)}{m_{+}I} \\ d_{+}power_{+}e\left( t \right) =& d_{+}R d_{+}power_{+}f\left( t \right) \\ s_{+}power_{+}e\left( t \right) =& \frac{s_{+}q\left( t \right)}{s_{+}C} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} =& s_{+}power_{+}f\left( t \right) \\ s3_{+}power_{+}e\left( t \right) =& \frac{\left( s3_{+}q\left( t \right) \right)^{3}}{s3_{+}C} \\ \frac{\mathrm{d} s3_{+}q\left( t \right)}{\mathrm{d}t} =& s3_{+}power_{+}f\left( t \right) \\ f_{+}power_{+}e\left( t \right) =& f_{+}F\left( t \right) \\ s_{+}power_{+}f\left( t \right) =& d_{+}power_{+}f\left( t \right) \\ s_{+}power_{+}f\left( t \right) =& m_{+}power_{+}f\left( t \right) \\ s_{+}power_{+}f\left( t \right) =& s3_{+}power_{+}f\left( t \right) \\ s_{+}power_{+}f\left( t \right) =& - f_{+}power_{+}f\left( t \right) \\ 0 =& - f_{+}power_{+}e\left( t \right) + d_{+}power_{+}e\left( t \right) + m_{+}power_{+}e\left( t \right) + s3_{+}power_{+}e\left( t \right) + s_{+}power_{+}e\left( t \right) \end{align} \]

The equations of the model above were not simplified and represent the bond-graphs equations directly derived from the connections. Then, we simplify the DAE above and obtain the ODE. The simplified set of equations is required for using solve from DifferentialEquations.jl

@named sys = simplifysys(mdl)

\[ \begin{align} \frac{\mathrm{d} m_{+}power_{+}f\left( t \right)}{\mathrm{d}t} =& \frac{\frac{ - \left( s3_{+}q\left( t \right) \right)^{3}}{s3_{+}C} + \frac{ - s_{+}q\left( t \right)}{s_{+}C} - d_{+}R m_{+}power_{+}f\left( t \right) + f_{+}F\left( t \right)}{m_{+}I} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} =& m_{+}power_{+}f\left( t \right) \\ \frac{\mathrm{d} s3_{+}q\left( t \right)}{\mathrm{d}t} =& m_{+}power_{+}f\left( t \right) \end{align} \]

Print system states

states(sys)

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 m₊power₊f(t)
 s₊q(t)
 s3₊q(t)

Print system parameters

parameters(sys)

5-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 d₊R
 f₊F(t)
 m₊I
 s3₊C
 s₊C

Print the system simplified equations

equations(sys)

By comparison, we can see that the equations above are the same as the canonical form of the harmonic oscillator presented at the beginning.

We can generate LaTeX code from the equations

latexify(equations(sys))

\begin{align} \frac{\mathrm{d} m{+}power{+}f\left( t \right)}{\mathrm{d}t} =& \frac{\frac{ - \left( s3{+}q\left( t \right) \right)^{3}}{s3{+}C} + \frac{ - s{+}q\left( t \right)}{s{+}C} - d{+}R m{+}power{+}f\left( t \right) + f{+}F\left( t \right)}{m{+}I} \ \frac{\mathrm{d} s{+}q\left( t \right)}{\mathrm{d}t} =& m{+}power{+}f\left( t \right) \ \frac{\mathrm{d} s3{+}q\left( t \right)}{\mathrm{d}t} =& m{+}power_{+}f\left( t \right) \end{align}

Simulate the system

Unforced case

Now, we are going to consider only the unforced case. Therefore, $F(t) = 0$. Thus, we update the force function on the system.

sys_unforced = substitute(sys, Dict(f.F => 0.0))

\[ \begin{align} \frac{\mathrm{d} m_{+}power_{+}f\left( t \right)}{\mathrm{d}t} =& \frac{\frac{ - \left( s3_{+}q\left( t \right) \right)^{3}}{s3_{+}C} + \frac{ - s_{+}q\left( t \right)}{s_{+}C} - d_{+}R m_{+}power_{+}f\left( t \right)}{m_{+}I} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} =& m_{+}power_{+}f\left( t \right) \\ \frac{\mathrm{d} s3_{+}q\left( t \right)}{\mathrm{d}t} =& m_{+}power_{+}f\left( t \right) \end{align} \]

Define the simulation time

tspan = (0, 10)

(0, 10)

When defining the ODEProblem as in the ModelingToolkit, we can define the initial value of the states.

prob = ODEProblem(sys_unforced, [m.power.f => 0.0, s.q => 1.0, s3.q => 1.0], tspan)

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 10)
u0: 3-element Vector{Float64}:
 0.0
 1.0
 1.0

Use the solve to simulate the system.

sol = solve(prob)
plot(sol, xlabel = "Time", ylabel = "Amplitude")

Forced case

Now, we are going to consider the forced case, where $F(t) = sin(t)$. Thus, we update the force function on the system.

sys_forced = substitute(sys, Dict(f.F => γ * cos(ω * t)))

\[ \begin{align} \frac{\mathrm{d} m_{+}power_{+}f\left( t \right)}{\mathrm{d}t} =& \frac{\frac{ - \left( s3_{+}q\left( t \right) \right)^{3}}{s3_{+}C} + \frac{ - s_{+}q\left( t \right)}{s_{+}C} + 8 \cos\left( 0.5 t \right) - d_{+}R m_{+}power_{+}f\left( t \right)}{m_{+}I} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} =& m_{+}power_{+}f\left( t \right) \\ \frac{\mathrm{d} s3_{+}q\left( t \right)}{\mathrm{d}t} =& m_{+}power_{+}f\left( t \right) \end{align} \]

Define the simulation time

tspan = (0, 10)

(0, 10)

When defining the ODEProblem as in the ModelingToolkit we can define the initial value of the states.

prob = ODEProblem(sys_forced, [m.power.f => 0.0, s.q => 0.0], tspan)

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 10)
u0: 3-element Vector{Float64}:
 0.0
 0.0
 0.0

Use the solve to simulate the system.

sol = solve(prob)
plot(sol, xlabel = "Time", ylabel = "Amplitude")

This page was generated using Literate.jl.