Reference

The main juliajim module

Useful routines for Harmonic Balance. Also includes Chebyshev Galerkin utilites.

ACT(yin, h::hTypes, N::Int64, dir::Symbol)

Routine for the alternating Chebyshev Transform. Uses AFT internally.

Arguments

• yin :
• h::hTypes :
• N::Int64 :
• dir::Symbol :

AFT(yin, h::hTypes, N::Int64, dir::Symbol);

Routine to conduct Time-to-Frequency and Frequency-to-Time transforms for the general Multi frequency case.

Arguments

• yin::VecOrMat{Float64}: (N^C, Ny) if time data, and (Nhc, Ny) if frequency data.
• h::hTypes: (Nh, C) list of harmonics.
• N::Int64: Number of time samples for AFT.
• dir::Symbol: :t2f for time-to-frequency, and :f2t for frequency-to-time.

Examples

• Single Frequency Case:

    N = 8;
    t = (0:N-1)*2π/N;
    yt = cos.(t) + 3sin.(2t) .+ 4;
    yt = [yt yt];

    h = collect(0:3)[:,:];
    Nhc = sum(all(h.==0, dims=2) + 2*any(h .!= 0, dims=2));

    yf = AFT(yt, h, N, :t2f);
    YT = AFT(yf, h, N, :f2t);

• 2-Frequency Case:

    t1 = repeat(t, 1, N);
    t2 = repeat(t', N, 1);

    yt = cos.(t1) + 3sin.(t2) .+ 4;
    h = [0 0;1 0;0 1;1 1];
    Nhc = sum(all(h.==0, dims=2) + 2*any(h .!= 0, dims=2));

    yf = AFT([yt[:] yt[:]], h, N, :t2f);
    YT = AFT(yf, h, N, :f2t);

DCHEB(h::hTypes)

Returns the Chebyshev Differentiation matrix.

Arguments

• h::hTypes :

DFOUR(h::hTypes, ws=nothing)

Returns the Fourier differentiation matrix (useful for HB representations).

Arguments

• h::hTypes : (H, C) Harmonic indices
• ws::Array{Float64} : (C, 1) (optional) list of frequencies

Outputs

• D : (Nhc, Nhc) Harmonic differentiation matrix
• dDdw : (Nhc, Nhc) Harmonic differentiation matrix (derivative wrt ws).

DFOUR!(h::hTypes, D, dDdw=nothing, ws=nothing)

Bang version of DFOUR.

Arguments

• h::hTypes : (H, C) Harmonic indices
• D : (Nhc, Nhc) Harmonic differentiation matrix
• dDdw : (Nhc, Nhc) Harmonic differentiation matrix (derivative wrt ws).
• ws::Array{Float64} : (C, 1) (optional) list of frequencies

FSEVAL(h::hTypes, t::Union{StepRangeLen{Float64}, VecOrMat{Float64}}, U::Union{Nothing,Matrix{Float64}}=nothing)

Description

Evaluates the Fourier series at required time points. If Fourier coefficients are not provided, it returns the mapping matrix.

Arguments

• h::hTypes :
• t::Union{StepRangeLen{Float64} :
• VecOrMat{Float64}} :
• U::Union{Nothing :
• Matrix{Float64}} : (default nothing)

Outputs

• ut : In case U is provided, this is outputted first.
• J : The mapping matrix such that J*U gives the desired temporal evaluates.

HARMONICSTIFFNESS!(E, dEdw, M, D, K, ws, h)

Bang version of HARMONICSTIFFNESS (preallocate E, dEdw).
Advantage is that nothing can be passed to avoid needless computations.

Arguments

• E: (NdNhc, NdNhc) Harmonic Stiffness
• dEdw: C-Vector of (NdNhc, NdNhc) Representing Frequency-gradients of E
• 'M, D, K': (Nd,Nd) Mass, Damping, Stiffness Matrices.
• 'ws::Array{Float64}': (C,1) list of frequencies
• 'h': (H, C) harmonic indices

Examples

• Single Frequency Case

M = I(2);
K = [2 -1;-1 2];
D = 0.001.*K + 0.01.*M;

h = [0;1;2;3];
Nhc = sum(all(h.==0, dims=2) + 2*any(h .!= 0, dims=2));
E = zeros(2*Nhc, 2*Nhc);
dEdw = zeros(2*Nhc, 2*Nhc);
ws = 1;
HARMONICSTIFFNESS!(E, dEdw, M, D, K, ws, h);

• 2-Frequency Case

Nhmax = 3;
ws = [1; π];
h = HSEL(Nhmax, ws);
Nhc = sum(all(h.==0, dims=2) + 2*any(h .!= 0, dims=2));
E = zeros(2*Nhc, 2*Nhc);
dEdw = zeros(2*Nhc, 2*Nhc);
HARMONICSTIFFNESS!(E, dEdw, M, D, K, ws, h);

HARMONICSTIFFNESS(M, D, K, ws::Array{Float64}, h)

Computes the harmonic stiffness for HB simulations.

Arguments

• 'M, D, K': (Nd,Nd) Mass, Damping, Stiffness Matrices.
• 'ws::Array{Float64}': (C,1) list of frequencies
• 'h': (H, C) harmonic indices

Outputs

• E: (NdNhc, NdNhc) Harmonic Stiffness
• dEdw: C-Vector of (NdNhc, NdNhc) Representing Frequency-gradients of E

Examples

• Single Frequency Case

M = I(2);
K = [2 -1;-1 2];
D = 0.001.*K + 0.01.*M;

h = [0;1;2;3];
ws = 1;
E, dEdw = HARMONICSTIFFNESS(M, D, K, ws, h);

• 2-Frequency Case

Nhmax = 3;
ws = [1; π];
h = HSEL(Nhmax, ws);
E, dEdw = HARMONICSTIFFNESS(M, D, K, ws, h);

HINDS(Ndofs::Int64, h::hTypes)

Arguments

• Ndofs::Int64 :
• h::hTypes :

Outputs

• zinds
• hinds
• rinds0
• rinds
• iinds

HSEL(Nhmax::Int64, ws::Any=1, hcr::Int=1)

Selection of Harmonic indices based on different criteria. First C+1 rows are [zeros(1,C); I(C)].

Implemented criteria are,

• hcr=1: Σᵢ |hᵢ| < Nhmax & ∑ᵢ h_i ≥ 0

Arguments

• Nhmax::Int64: Max harmonic for truncation.
• ws=1: [opt] List of frequencies
• hcr::Int=1: [opt] Truncation criterion. See above.

Examples

C = 2;
Nhmax = 3;

h = HSEL(Nhmax, 1:C)

NHC(h::hTypes)

Returns the number of harmonics in a chosen harmonic set h.

Arguments

• h::hTypes :

Description

Arguments

• U :
• h :
• Hmax : (default nothing)
• D : (default nonothing)
• L : (default nonothing)

PRODMAT_FOUR(U::VecOrMat{Float64}, h::hTypes, Hmax=nothing, D=nothing, L=nothing)

Description

Returns the Fourier product matrix. Currently only implemented for C=1.

Arguments

• U::VecOrMat{Float64} :
• h::hTypes :
• Hmax : (default nothing)
• D : (default nothing)
• Lb : (default nothing)

Useful types and routines for numerical continuation.

myNLSoln(up::nvTypes=nothing; J::nmTypes=nothing, Jp::nvTypes=nothing)

Constructor for myNLSoln. Can specify one or all the arguments (point, jacobian, paramjac).

Arguments

• up::nvTypes : (default nothing) Solution point
• J::nmTypes : (default nothing) Jacobian
• Jp::nvTypes : (default nothing) paramjac
• save_jacs::Bool : (default false) Whether or not to store the Jacobians.

If both J and Jp are provided, the unit tangent dupds is computed using nullspace([J Jp])[:,1].

myNLSoln

Struct storing the solution, Jacobian, and (unit) tangent.

Fields

• up::nvTypes: Solution point ([u;p])
• J::nmTypes: Jacobian dr/du
• Jp::nvTypes: Parameter Jacobian dr/dp
• dupds::nvTypes: Tangent vector (has to be vector only)

Base.:-(v1::myNLSoln, v2::myNLSoln)

Overloading the subtraction operator to take the difference between two myNLSoln

Base.getproperty(sols::Vector{myNLSoln}, sym::Symbol)

This will allow accessing elements of myNLSoln from vectors of structs.

Arguments

• sols::Vector{myNLSoln} :
• sym::Symbol :

Base.getproperty(sol::myNLSoln, sym::Symbol)

This will allow accessing the solution (u) and parameter (p) from myNLSoln easily

Arguments

• sol::myNLSoln :
• sym::Symbol :

Base.show(io::IO, p::myNLSoln)

Overloading base show function to display myNLSoln object.

CONTINUATE(u0, fun, ps, dp; pkwargs...)

Continuation routine. Solves the bordered problem with residue drawn from EXTRESFUN!.

Arguments

• u0::Vector{Float64} : Initial guess for first point
• fun : NonlinearFunction object. Will call:
  • fun.f(du, u, p) for residue r. p must be scalar.
  • fun.jac(J, u, p) for jacobian dr/du. Uses ForwardDiff.jacobian! if fun.jac is nothing.
  • fun.paramjac(Jp, u, p) for jacobian dr/dp. Uses ForwardDiff.jacobian! if fun.paramjac is nothing.
• ps::Vector{Float64} : Specify range of p for continuation.
• dp::Float64 : Specify first step in Δp units. (will be rescaled if necessary)

Optional Arguments

• parm::Symbol : (default :arclength) Arclength parameterization. See EXTRESFUN!.
• nmax::Int64 : (default 1000) Maximum number of steps.
• dpbnds::Union{Nothing,Vector{Float64}} : (default [dp/5,5dp]) Bounds for the step length dp. (Rescaled as appropriate)
• save_jacs::Bool : (default false) Specify whether or not to save the Jacobians in the output. Only solution and (unit) tangent are saved if false.
• verbosity::Int : (default 1) Verbosity levels:
  • 0: Suppress all messages.
  • 1: Display Steps, Continuation Progress.
  • 2: Display Iteration Information for each step also.

For Scaling of Unknowns (recommended to get evenly spaced points on response curve)

• Dsc::Union{Symbol,Nothing,Vector{Float64}} : (default :none) "Dscaling" used to scale the unknowns. The arc length constraint is applied in the scaled space (uₛ=uₚₕ./Dsc).
  • If set to :auto, it uses the absolute of the first converged solution as the initial Dsc vector. Zero entries are replaced with minDsc.
  • If set to :none, it fixes Dsc to a vector of ones and doesn't dynamically adapt it. (this forces DynScale to false).
  • If set to :ones, it fixes Dsc to a vector of ones but dynamically continues to scale.
• DynScale::Bool : (default true) Whether or not to dynamically adapt the Dsc vector. Each entry is allowed to grow or shrink by a maximum factor of 2 in each step if true.
• minDsc::Float64 : (default eps()^(4//5)=3e-13) Minimum value for Dscale.

For Step Length Adaptation. Currently set as dsₙ = dsₒ * xi, where xi = clamp((itopt/itns)^nxi, xirange[1], xirange[2]).

• itopt::Union{Symbol,Int} : (default :auto) Optimal number of
• nxi::Float64 : (default 0.0) Switches off adaptation by default.
• xirange::Vector{Float64} : (default [0.5,2.])

Returns

• vector{myNLsoln} representing solutions.
• vector{Int} representing iterations taken.
• vector{Float64} representing step sizes.
• vector{Float64} representing adaptation parameter ξ.
• vector{Float64} representing the final Dscale matrix.

EXTRESFUN!(up::nvTypes, fun, sol0::myNLSoln, ds::Float64; parm::Symbol=:riks, Dsc::Vector{Float64}, dup::nvTypes=nothing, Jf::nmTypes=nothing)

Returns the bordered/extended residue function. Appends residue with required arclength constraint.

Note

The function involves very naive autodiff calls. Does not respect or try to detect jacobian sparsity. For now it is best to provide analytical jacobians for very large problems.

Arguments

• up::nvTypes: Solution point for evaluation

• fun: Nonlinear function specified in a fashion that can be called as

  • fun.f(du,up[1:end-1],up[end]) for the residue.

  • fun.jac(J,up[1:end-1],up[end]) for the jacobian (uses ForwardDiff.jacobian! if fun.jac is nothing)

  • fun.paramjac(Jp,up[1:end-1],up[end]) for the parameter jacobian (uses ForwardDiff.jacobian! if fun.paramjac is nothing)

  • Development Done based on NonlinearFunction from NonlinearSolve

• sol0::myNLSoln: Previous solution point.

• ds::Float64: Step size.

• parm::Symbol=:riks: Arclength parameterization. Defaults to :riks. Possible are:

  • :riks: Riks' or normal parameterization.

  • :arclength: Arclength parameterization.

• Dsc::Vector{Float64}: Scaling vector. Recommended to be of the order of magnitude of the expected unknowns. Size same as up.

• dup::nvTypes=nothing: Residue vector (for iip evaluation)

• Jf::nmTypes=nothing: Jacobian matrix (for iip evaluation)

Returns

• dup::nvTypes=nothing: Residue vector (for iip evaluation)

• Jf::nmTypes=nothing: Jacobian matrix (for iip evaluation)

Useful types and routines for dynamic analyses of MDOF problems.

MDOFGEN

Struct storing the generic MDOF system.

Fields

• M::Matrix{Float64}: Mass matrix
• C::Matrix{Float64}: Damping Matrix
• K::Matrix{Float64}: Stiffness Matrix
• L: Displacement Transform Matrix (can be empty)
• Ndofs::Int64: Number of DOFs
• NLTs::Vector{NONLINEARITY}: Vector of nonlinearities present (see NONLINEARITY)

NONLINEARITY

Struct storing a nonlinearity.

Fields

• type::Symbol: Type of nonlinearity. One of:

  1. :Inst (for instantaneous nonlinearity)

  2. :Hyst (for hysteretic nonlinearity)

  3. :Fdom (for nonlinearity defined in frequency domain)

     Out of these, only the first 2 can be used for transient simulations.

• func::Function: Function handle for evaluating nonlinearity. Has signature

  • (t, u, udot) for :Inst
  • (t, u, up, sp) for :Hyst, where s is some internal state (sp is s at t-Δt)

• L::Matrix{Float64}: Selection Matrix

• Lf: Force shape matrix (can be empty)

Description

Add a nonlinearity to the MDOFGEN object.

Arguments

• m::MDOFGEN : MDOFGEN object
• type::Symbol : One of :Inst (Instantaneous nonlinearity), :Hyst (Hysteretic nonlinearity), :Fdom (Frequency Domain nonlinearity), :Dla (Dynamic Lagrangian)
• func::Function : Nonlinearity function returning the nonlinear force. Different signatures for different types:
  • :Inst : (t,u,ud) returning (f, dfdu, dfdud)
  • :Hyst : (t,u,uprev,fprev) returning (f, dfdu, dfduprev, dfdfprev)
  • :Fdom : (Unl, h, N) returning (F, dFdU, dFdw)
  • :Dlag : Not implemented yet
• L :
• Lf : (default nothing)

Description

Arguments

• U :
• U0 : Vector{Float64} or Vector{Vector{Float64}}
• resfun! : function!(u, Res, Jac, Jacp)
• R : (default nothing)
• J : (default nothing)
• Jp : (default nothing)

Description

Arguments

• Uwxa :
• m::MDOFGEN :
• Fl : (Nd*Nhc) Forcing vector. Zero harmonics are used as static loads. Harmonic portions are used to provide phase constraint.
• h :
• N::Int64 :
• tol::Float64 : (default eps()^(4//5))
• atype::Symbol : (default :H1)
• ashape::Union{Int64 :
• Vector{Float64}} : (default 1)
• R : (default nothing)
• dRdUwx : (default nothing)
• dRda : (default nothing)

Description

Residue function that can be used for Harmonic Balance forced response analysis.

Arguments

• Uw :
• m :
• Fl :
• h :
• N :
• tol :
• dRdU :
• dRdw :

Description

HB residue with amplitude constraint.

Arguments

• Ufw :
• m::MDOFGEN :
• A::Float64 : Amplitude level (also see atype below)
• Fl :
• h :
• N::Int64 :
• atype::Symbol : (default :H1) One of :H1, :RMS. Type of amplitude measure.
• ashape : (default 1) Shape to dof for applying amplitude. If Integer, interpreted as Dof. If vector, taken as shape.
• tol::Float64 : (default eps()^(4//5))
• R : (default nothing)
• dRdUf : (default nothing)
• dRdw : (default nothing)

Description

Arguments

• m::MDOFGEN :
• T0 :
• T1 :
• dt :
• U0 :
• Ud0 :
• Fex :

Description

Evaluate the nonlinearities (and Jacobian) for one period and return the harmonics

Arguments

• Uw :
• m::MDOFGEN :
• h :
• N::Int64 : [optional]
• tol :
• FNL :
• dFNLdU :
• dFNLdw :

Description

Evaluate the nonlinearities in the time domain for a single time instant.

Arguments

• m::MDOFGEN : MDOFGEN object.
• t : Time
• U : (m.Ndofs) Displacements
• Ud : (m.Ndofs) Velocities
• tp : (default nothing) Previous instant time
• Up : (default nothing) (m.Ndofs) Displacements at previous instant.
• Udp : (default nothing) (m.Ndofs) Velocities at previous instant.
• Sp : (default nothing) (length(m.NLTs)) Vector of vectors of internal states at previous instant. Each nonlinearity can have its own set of internal states.

Returns

• F : (m.Ndofs) force vector
• dFdU : (m.Ndofs, m.Ndofs) force-displacement derivative matrix
• dFdUd : (m.Ndofs, m.Ndofs) force-velocity derivative matrix
• S : (length(m.NLTs)) vector of vectors of internal states

Returns the quasi periodic forced response residue. Also allows for specifying a constrained form. If provided, the matrix cL and vector cLb are interpreted such that the unknown vector gets written as $U = {}^cL Û + {}^cL_b,$ where Û is of size (D Nhc+C,1), containing the frequencies also in the end.

The matrix {}^cL is used to project the residue vector {}^cL^T R. By choosing its size appropriately, one can make the system square.

The last element in U is always interpreted as the continuation parameter.

Arguments

• Uw :
• m::MDOFGEN :
• Fl::Union{Vector :
• Matrix :
• SparseMatrixCSC :
• Nothing} :
• h :
• N::Int64 :
• tol::Float64 : (default eps()^(4//5))
• R : (default nothing)
• dRdU : (default nothing)
• dRdw : (default nothing)
• cL : (default nothing)
• cLb : (default nothing)
• U0 : (default nothing) Solution to deflate from. 2-norm is used for deflation.

Some routines for Nonlinear Dynamics.

Description

NORMALFORMFIT fits a quadratic normal form to the given residue function (in R2) numerically. It starts with a cloud of randomly distributed points with given amplitude standard deviation and fits a general cubic model such that the dynamics may be approximated as: xdot = A x + B x x + C x x x . Then the quadratic terms are eliminated by posing a near identity transformation from x to z as x = z + H z z . The simplified normal form equations in terms of z are zdot = A z + C3 z z z .

A good reference for this is sec. 19.10 in Wiggins (1990).

Arguments

• xyfun : (xy) Taking in a 2-vector and returning a 2-residue vector.

	Must have equilibrium at origin.

• damp : If scalar, standard deviation for amplitude. If vector, [dmin, d_max] for the amplitude.
• Npts : Number of points to be used for fit.
• constrained : (default true)

Outputs

• A : (2,2) The linear Jacobian matrix.
• B2n : (2,3) The quadratic coefficient matrix such that the nonlinear influence is written as B2[z1^2; z1z2; z2^2]
• C3n : (2,4) The cubic coefficient matrix such that the nonlinear influence is written as C3[z1^3; z1^2z2; z1*z2^2; z2^3].
• Hn : (2,3) The cubic near identity transformation matrix such that the near identity transformation can be written as x = z + Hn [z1^2; z1*z2; z2^2]
• C3 : (2,2,2,2) The cubic coefficient matrix in tensor form.
• H : (2,2,2) The quadratic near identity transformation coefficients in tensor form.