System functions¶

The group of system functions comprises functions for the setting up, solving, and elimination of systems of equations. The functions are separated in two groups:

Static system functions
Dynamic system functions

Static system functions concern the linear system of equations

\[\mathbf{K} \mathbf{a} = \mathbf{f}\]

where \(\mathbf{K}\) is the global stiffness matrix and \(\mathbf{f}\) is the global load vector. Often used static system functions are assem and solveq. The function assem assembles the global stiffness matrix and solveq computes the global displacement vector \(\mathbf{a}\) considering the boundary conditions. It should be noted that \(\mathbf{K}\), \(\mathbf{f}\), and \(\mathbf{a}\) also represent analogous quantities in systems other than structural mechanical systems. For example, in a heat flow problem \(\mathbf{K}\) represents the conductivity matrix, \(\mathbf{f}\) the heat flow, and \(\mathbf{a}\) the temperature.

Dynamic system functions are related to different aspects of linear dynamic systems of coupled ordinary differential equations according to

\[\mathbf{C} \dot{\mathbf{a}} + \mathbf{K} \mathbf{a} = \mathbf{f}\]

for first-order systems and

\[\mathbf{M} \ddot{\mathbf{a}} + \mathbf{C} \dot{\mathbf{a}} + \mathbf{K} \mathbf{a} = \mathbf{f}\]

for second-order systems. First-order systems occur typically in transient heat conduction and second-order systems occur in structural dynamics.

Static system functions¶

The group of static system functions comprises functions for setting up and solving the global system of equations. It also contains a function for eigenvalue analysis, a function for static condensation, a function for extraction of element displacements from the global displacement vector and a function for extraction of element coordinates.

assem¶

Purpose:: Assemble element matrices.
Syntax:

K = assem(edof, K, Ke)
[K, f] = assem(edof, K, Ke, f, fe)

Description:

assem adds the element stiffness matrix \(\mathbf{K}^e\), stored in Ke, to the structure stiffness matrix \(\mathbf{K}\), stored in K, according to the topology matrix edof.

The element topology matrix edof is defined as

edof\(= [el \quad \underbrace{dof_1\quad dof_2\quad\ldots \quad dof_{ned}}_{\text{global dof.}} ]\)

where the first column contains the element number, and columns 2 to \((ned+1)\) contain the corresponding global degrees of freedom (\(ned\) = number of element degrees of freedom).

In the case where the matrix \(\mathbf{K}^e\) is identical for several elements, assembling of these can be carried out simultaneously. Each row in edof then represents one element, i.e. \(nel\) is the total number of considered elements.

Edof\(= \left. \left[ \begin{array}{c} el_1 \\ el_2 \\ \vdots \\ el_{nel} \end{array} \quad \begin{array}{cccccc} dof_1 & dof_2 & \cdots & \cdots & \cdots & dof_{ned} \\ dof_1 & dof_2 & \cdots & \cdots & \cdots & dof_{ned} \\ \vdots & \vdots & & & & \vdots \\ dof_1 & dof_2 & \cdots & \cdots & \cdots & dof_{ned} \end{array} \right] \right\} \text{one row for each element}\)

If \(\mathbf{fe}\) and \(\mathbf{f}\) are given in the function, the element load vector \(\mathbf{f}^e\) is also added to the global load vector \(\mathbf{f}\).

coordxtr¶

Purpose:: Extract element coordinates from a global coordinate matrix.
Syntax:

[Ex, Ey, Ez] = coordxtr(Edof, Coord, Dof, nen)

Description:

coordxtr extracts element nodal coordinates from the global coordinate matrix Coord for elements with equal numbers of element nodes and dof’s.

Input variables are the element topology matrix Edof, defined in assem, the global coordinate matrix Coord, the global topology matrix Dof, and the number of element nodes \(nen\) in each element.

Coord\(= \begin{bmatrix} x_1 & y_1 & [z_1] \\ x_2 & y_2 & [z_2] \\ x_3 & y_3 & [z_3] \\ \vdots & \vdots & \vdots \\ x_n & y_n & [z_n] \end{bmatrix}\) \(\qquad\) Dof\(= \begin{bmatrix} k_1 & l_1 & \ldots & m_1 \\ k_2 & l_2 & \ldots & m_2 \\ k_3 & l_3 & \ldots & m_3 \\ \vdots & \vdots & \ldots & \vdots \\ k_n & l_n & \ldots & m_n \end{bmatrix}\) \(\qquad\) nen\(= [\;nen\;]\)

The nodal coordinates defined in row \(i\) of Coord correspond to the degrees of freedom of row \(i\) in Dof. The components \(k_i\), \(l_i\) and \(m_i\) define the degrees of freedom of node \(i\), and \(n\) is the number of global nodes for the considered part of the FE-model.

The output variables Ex, Ey, and Ez are matrices defined according to

Ex\(= \begin{bmatrix} x_1^1 & x_2^1 & x_3^1 & \ldots & x_{nen}^1 \\ x_1^2 & x_2^2 & x_3^2 & \ldots & x_{nen}^2 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_1^{nel} & x_2^{nel} & x_3^{nel} & \ldots & x_{nen}^{nel} \end{bmatrix}\)

where row \(i\) gives the \(x\)-coordinates of the element defined in row \(i\) of Edof, and where \(nel\) is the number of considered elements.

The element coordinate data extracted by the function coordxtr can be used for plotting purposes and to create input data for the element stiffness functions.

eigen¶

Purpose:: Solve the generalized eigenvalue problem.
Syntax:

L = eigen(K, M)
L = eigen(K, M, b)
[L, X] = eigen(K, M)
[L, X] = eigen(K, M, b)

Description:

eigen solves the eigenvalue problem

\[| \mathbf{K} - \lambda \mathbf{M} | = 0\]

where \(\mathbf{K}\) and \(\mathbf{M}\) are square matrices. The eigenvalues \(\lambda\) are stored in the vector \(\mathbf{L}\) and the corresponding eigenvectors in the matrix \(\mathbf{X}\).

If certain rows and columns in matrices \(\mathbf{K}\) and \(\mathbf{M}\) are to be eliminated in computing the eigenvalues, b must be given in the function. The rows (and columns) that are to be eliminated are described in the vector b defined as

b\(=\begin{bmatrix} dof_1 \\ dof_2 \\ \vdots \\ dof_{nb} \end{bmatrix}\)

The computed eigenvalues are given in order ranging from the smallest to the largest. The eigenvectors are normalized so that

\[\mathbf{X}^T \mathbf{M} \mathbf{X} = \mathbf{I}\]

where \(\mathbf{I}\) is the identity matrix.

extract_ed¶

Purpose:: Extract element nodal quantities from a global solution vector.
Syntax:

ed = extract_ed(edof, a)

Description:

The extract_ed function extracts element displacements or corresponding quantities \(\mathbf{a}^e\) from the global solution vector \(\mathbf{a}\), stored in a.

Input variables are the element topology matrix edof, defined in assem, and the global solution vector a.

The output variable

ed\(= (\mathbf{a}^e)^T\)

contains the element displacement vector.

If Edof contains more than one element, Ed will be a matrix

Ed\(= \begin{bmatrix} (\mathbf{a}^e)_1^T \\ (\mathbf{a}^e)_2^T \\ \vdots \\ (\mathbf{a}^e)_{nel}^T \end{bmatrix}\)

where row \(i\) gives the element displacements for the element defined in row \(i\) of Edof, and \(nel\) is the total number of considered elements.

Example:

For the two-dimensional beam element, the extract_ed function will extract six nodal displacements for each element given in Edof, and create a matrix Ed of size \((nel × 6)\).

Ed\(= \begin{bmatrix} u_1 & u_2 & u_3 & u_4 & u_5 & u_6 \\ u_1 & u_2 & u_3 & u_4 & u_5 & u_6 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ u_1 & u_2 & u_3 & u_4 & u_5 & u_6 \end{bmatrix}\)

red¶

Purpose:: Reduce the size of a square matrix by omitting rows and columns.
Syntax:

B = red(A, b)
[B, b] = red(A, b)

Description:

B = red(A, b) reduces the square matrix A to a smaller matrix B by omitting rows and columns of A. The indices for rows and columns to be omitted are specified by the column vector b.

Example:

Assume that the matrix A is defined as

A\(= \begin{bmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \end{bmatrix}\)

and b as

b\(= \begin{bmatrix} 2 \\ 4 \end{bmatrix}\)

The statement B = red(A, b) results in the matrix

B\(= \begin{bmatrix} 1 & 3 \\ 9 & 11 \end{bmatrix}\)

solveq¶

Purpose:: Solve equation system.
Syntax:

a = solveq(K, f)
a = solveq(K, f, bcP, bc)
[a, r] = solveq(K, f, bc)

Description:

The function solveq solves the equation system

\[\mathbf{K}\;\mathbf{a} = \mathbf{f}\]

where \(\mathbf{K}\) is a matrix and \(\mathbf{a}\) and \(\mathbf{f}\) are vectors.

The matrix \(\mathbf{K}\) and the vector \(\mathbf{f}\) must be predefined. The solution of the system of equations is stored in a vector \(\mathbf{a}\) which is created by the function.

If some values of \(\mathbf{a}\) are to be prescribed, the row number and the corresponding values are given in the boundary condition matrix

bc\(= \left[ \begin{array}{c} dof_1 \\ dof_2 \\ \vdots \\ dof_{nbc} \end{array} \quad \begin{array}{c} u_1 \\ u_2 \\ \vdots \\ u_{nbc} \end{array} \right]\)

where the first column contains the row numbers and the second column the corresponding prescribed values.

If r is given in the function, support forces are computed according to

mathbf{r}\(= \mathbf{K}\;\mathbf{a} - \mathbf{f}\)

statcon¶

Purpose:: Reduce system of equations by static condensation.
Syntax:

[K1, f1] = statcon(K, f, b)

Description:

statcon reduces a system of equations

\[\mathbf{K}\;\mathbf{a} = \mathbf{f}\]

by static condensation.

The degrees of freedom to be eliminated are supplied to the function by the vector

b\(= \begin{bmatrix} dof_1 \\ dof_2 \\ \vdots \\ dof_{nb} \end{bmatrix}\)

where each row in b contains one degree of freedom to be eliminated.

The elimination gives the reduced system of equations

\[\mathbf{K}_1\;\mathbf{a}_1 = \mathbf{f}_1\]

where \(\mathbf{K}_1\) and \(\mathbf{f}_1\) are stored in K1 and f1 respectively.

Dynamic system functions¶

The group of system functions comprises functions for solving linear dynamic systems by time stepping or modal analysis, functions for frequency domain analysis, etc.

dyna2¶

Compute the dynamic solution to a set of uncoupled second-order differential equations.

Syntax:

X = dyna2(w2, xi, f, g, dt)

Description:

dyna2 computes the solution to the set

\[\ddot{x}_i + 2 \xi_i \omega_i \dot{x}_i + \omega^{2}_i x_i = f_i g(t), \qquad i=1,\ldots,m\]

of differential equations, where \(g(t)\) is a piecewise linear time function.

The vectors w2, xi, and f contain the squared circular frequencies \(\omega_i^2\), the damping ratios \(\xi_i\), and the applied forces \(f_i\), respectively.

The vector g defines the load function in terms of straight line segments between equally spaced points in time. This function may have been formed by the command gfunc.

The dynamic solution is computed at equal time increments defined by dt. Including the initial zero vector as the first column vector, the result is stored in the \(m \times n\) matrix X, where \(n-1\) is the number of time steps.

Note

The accuracy of the solution is not a function of the output time increment dt, since the command produces the exact solution for straight line segments in the loading time function.

See also:: gfunc

dyna2f¶

Purpose:

Compute the dynamic solution to a set of uncoupled second-order differential equations.

Syntax:

Y = dyna2f(w2, xi, f, p, dt)

Description:

dyna2f computes the solution to the set of differential equations:

\[\ddot{x}_i + 2 \xi_i\omega_i \dot{x}_i + \omega^{2}_i x_i = f_i g(t), \qquad i=1,\ldots,m\]

The vectors w2, xi and f are the squared circular frequencies \(\omega_i^2\), the damping ratios \(\xi_i\), and the applied forces \(f_i\), respectively. The force vector p contains the Fourier coefficients \(p(k)\) formed by the command fft.

The solution in the frequency domain is computed at equal time increments defined by dt. The result is stored in the \(m \times n\) matrix Y, where m is the number of equations and n is the number of frequencies resulting from the fft command. The dynamic solution in the time domain is achieved by the use of the command ifft.

Example:

The dynamic solution to a set of uncoupled second-order differential equations can be computed by the following sequence of commands:

>> g = gfunc(G, dt);
>> p = fft(g);
>> Y = dyna2f(w2, xi, f, p, dt);
>> X = (real(ifft(Y.')))';

where it is assumed that the input variables G, dt, w2, xi and f are properly defined. Note that the ifft command operates on column vectors if Y is a matrix; therefore use the transpose of Y. The output from the ifft command is complex. Therefore use Y.' to transpose rows and columns in Y in order to avoid the complex conjugate transpose of Y.

The time response is represented by the real part of the output from the ifft command. If the transpose is used and the result is stored in a matrix X, each row will represent the time response for each equation as the output from the command dyna2.

See also:

gfunc, fft, ifft

fft¶

Purpose:

Transform functions in time domain to frequency domain.

Syntax:

p = fft(g)
p = fft(g, N)

Description:

The fft function transforms a time dependent function to the frequency domain.

The function to be transformed is stored in the vector g. Each row in g contains the value of the function at equal time intervals. The function represents a span \(-\infty \leq t \leq +\infty\); however, only the values within a typical period are specified by g.

The fft command can be used with one or two input arguments. If N is not specified, the number of frequencies used in the transformation is equal to the number of points in the time domain (i.e., the length of the variable g), and the output will be a vector of the same size containing complex values representing the frequency content of the input signal.

The scalar variable N can be used to specify the number of frequencies used in the Fourier transform. The size of the output vector in this case will be equal to N. It should be remembered that the highest harmonic component in the time signal that can be identified by the Fourier transform corresponds to half the sampling frequency. The sampling frequency is equal to \(1/dt\), where \(dt\) is the time increment of the time signal.

The complex Fourier coefficients \(p(k)\) are stored in the vector p, and are computed according to

\[p(k) = \sum_{j=1}^N x(j) \omega_N^{(j-1)(k-1)},\]

where

\[\omega_N = e^{-2 \pi i / N}.\]

Note

This is a MATLAB built-in function.

freqresp¶

Purpose:

Compute frequency response of a known discrete time response.

Syntax:

[Freq, Resp] = freqresp(D, dt)

Description:

freqresp computes the frequency response of a discrete dynamic system.

D is the time history function.
dt is the sampling time increment, i.e., the time increment used in the time integration procedure.
Resp contains the computed response as a function of frequency.
Freq contains the corresponding frequencies.

Example:

The result can be visualized by:

plot(Freq, Resp)
xlabel('frequency (Hz)')

or:

semilogy(Freq, Resp)
xlabel('frequency (Hz)')

The dimension of Resp is the same as that of the original time history function.

Note

The time history function of a discrete system computed by direct integration often behaves in an unstructured manner. The reason for this is that the time history is a mixture of several participating eigenmodes at different eigenfrequencies. By using a Fourier transform, however, the response as a function of frequency can be computed efficiently. In particular, it is possible to identify the participating frequencies.

gfunc¶

Purpose:

Form vector with function values at equally spaced points by linear interpolation.

Syntax:

[t, g] = gfunc(G, dt)

Description:

gfunc uses linear interpolation to compute values at equally spaced points for a discrete function \(g\) given by

\[\begin{split}G = \left[ \begin{array}{cc} t_1 & g(t_1)\\ t_2 & g(t_2)\\ \vdots \\ t_N & g(t_N) \end{array} \right],\end{split}\]

as shown in the figure above.

Function values are computed in the range \(t_1 \leq t \leq t_N\), at equal increments, dt being defined by the variable dt. The number of linear segments (steps) is \((t_N-t_1)/dt\). The corresponding vector t is also computed. The result can be plotted by using the command plot(t, g).

ifft¶

Purpose:

Transform function in frequency domain to time domain.

Syntax:

x = ifft(y)
x = ifft(y, N)

Description:

ifft transforms a function in the frequency domain to a function in the time domain.

The function to be transformed is given in the vector y. Each row in y contains Fourier terms in the interval \(-\infty \leq \omega \leq +\infty\).

The ifft command can be used with one or two input arguments. The scalar variable N can be used to specify the number of frequencies used in the Fourier transform. The size of the output vector in this case will be equal to N. See also the description of the command fft.

The inverse Fourier coefficients \(x(j)\), stored in the variable x, are computed according to

\[x(j) = \frac{1}{N} \sum_{k=1}^N y(k) \omega_N^{-(j-1)(k-1)},\]

where

\[\omega_N = e^{-2 \pi i / N}.\]

Note

This is a MATLAB built-in function.

See also:

fft

ritz¶

Purpose:

Compute approximative eigenvalues and eigenvectors by the Lanczos method.

Syntax:

L = ritz(K, M, f, m)
L = ritz(K, M, f, m, b)
[L, X] = ritz(K, M, f, m)
[L, X] = ritz(K, M, f, m, b)

Description:

ritz computes, by the use of the Lanczos algorithm, m approximative eigenvalues and m corresponding eigenvectors for a given pair of n-by-n matrices K and M and a given non-zero starting vector f.

If certain rows and columns in matrices \(\mathbf{K}\) and \(\mathbf{M}\) are to be eliminated in computing the eigenvalues, \(\mathbf{b}\) must be given in the command. The rows (and columns) to be eliminated are described in the vector \(\mathbf{b}\) defined as

\[\begin{split}\mathbf{b} = \begin{bmatrix} dof_1 \\ dof_2 \\ \vdots \\ dof_{nb} \end{bmatrix}\end{split}\]

Note

If the number of vectors, m, is chosen less than the total number of degrees-of-freedom, \(n\), only about the first m/2 Ritz vectors are good approximations of the true eigenvectors. Recall that the Ritz vectors satisfy the M-orthonormality condition

\[\mathbf{X}^T \mathbf{M} \mathbf{X} = \mathbf{I}\]

where \(\mathbf{I}\) is the identity matrix.

spectra¶

Purpose:

Compute seismic response spectra for elastic design.

Syntax:

s = spectra(a, xi, dt, f)

Description:

The spectra function computes the seismic response spectrum for a known acceleration history function.

The computation is based on the vector a, which contains an acceleration time history function defined at equal time steps.
The time step is specified by the variable dt.
The value of the damping ratio is given by the variable xi.
Output from the computation, stored in the vector s, is achieved at frequencies specified by the column vector f.

Example:

The following procedure can be used to produce a seismic response spectrum for a damping ratio \(\xi = 0.05\), defined at 34 logarithmically spaced frequency points. The acceleration time history a has been sampled at a frequency of 50 Hz, corresponding to a time increment dt = 0.02 between collected points:

freq = logspace(0, log10(2^(33/6)), 34);
xi = 0.05;
dt = 0.02;
s = spectra(a, xi, dt, freq');

The resulting spectrum can be plotted by the command:

loglog(freq, s, '*')

step1¶

Purpose:

Compute the dynamic solution to a set of first order differential equations.

Syntax:

[a,da] = step1(K, C, f, a0, bc, ip)
[a,da] = step1(K, C, f, a0, bc, ip, times)
[a,da, ahist, dahist] = step1(K, C, f, a0, bc, ip, times, dofs)

Description:

step1 computes at equal time steps the solution to a set of first order differential equations of the form:

\[\begin{split}\mathbf{C} \dot{\mathbf{a}} + \mathbf{K}\mathbf{a} = \mathbf{f}(x, t), \\ \mathbf{a}(0) = \mathbf{a}_0\end{split}\]

The command solves transient field problems. In the case of heat conduction, K and C represent the \(n \times n\) conductivity and capacity matrices, respectively. a is the temperature and da (i.e., \(\dot{\mathbf{a}}\)) is the time derivative of the temperature.

The matrix f contains the time-dependent load vectors. If no external loads are active, use [] for f. The matrix f is organized as follows:

f = [
time history of the load at dof_1
time history of the load at dof_2
...
time history of the load at dof_n
]

The dimension of f is:

(number of degrees-of-freedom) × (number of timesteps + 1)

The initial conditions are given by the vector a0 containing initial values of a.

The matrix bc contains the time-dependent prescribed values of the field variable a. If no field variables are prescribed, use [] for bc. The matrix bc is organized as follows:

bc = [
dof_1   time history of the field variable
dof_2   time history of the field variable
...
dof_m2  time history of the field variable
]

The dimension of bc is:

(number of dofs with prescribed field values) × (number of timesteps + 2)

The time integration procedure is governed by the parameters given in the vector ip defined as:

ip = [dt, T, alpha]

where dt specifies the length of the time increment, T is the total time, and alpha is the time integration constant. Frequently used values of alpha are:

alpha	Method
0	Forward difference; forward Euler
0.5	Trapezoidal rule; Crank-Nicholson
1	Backward difference; backward Euler

The computed values of a and da are stored in a and da, respectively. The first column contains the initial values, and the following columns contain the values for each time step. The dimension is:

(number of degrees-of-freedom) × (number of time steps + 1)

If the values are to be stored only for specific times, the parameter times specifies at which times the solution will be stored. The values are stored in a and da, one column for each requested time according to times. The dimension is then:

(number of degrees-of-freedom) × (number of requested times + 1)

If the history is to be stored in ahist and dahist for some degrees of freedom, the parameter dofs specifies for which degrees of freedom the history is to be stored. The computed time histories are stored in ahist and dahist, respectively, with one row for each requested degree of freedom. The dimension is:

(number of specified degrees of freedom) × (number of timesteps + 1)

The time history functions can be generated using the command gfunc. If all the values of the time histories of f or bc are kept constant, these values need to be stated only once. In this case, the number of columns in f is one and in bc two.

In most cases, only a few degrees-of-freedom are affected by the exterior load, and hence the matrix contains only a few non-zero entries. In such cases, it is possible to save space by defining f as sparse (a MATLAB built-in function).

Note

Reference: Bathe, K.J.: Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, New Jersey, pp. 511-514, 1982.

step2¶

Purpose:

Compute the dynamic solution to a set of second order differential equations.

Syntax:

[a, da, d2a] = step2(K, C, M, f, a0, da0, bc, ip)
[a, da, d2a] = step2(K, C, M, f, a0, da0, bc, ip, times)
[a, da, d2a, ahist, dahist, d2ahist] = step2(K, C, M, f, a0, da0, bc, ip, times, dofs)

Description:

step2 computes at equal time steps the solution to a set of second order differential equations of the form:

\[\begin{split}\mathbf{M} \ddot{\mathbf{a}} + \mathbf{C} \dot{\mathbf{a}} + \mathbf{K} \mathbf{a} = \mathbf{f}(x, t), \\ \mathbf{a}(0) = \mathbf{a}_0, \\ \dot{\mathbf{a}}(0) = \mathbf{v}_0.\end{split}\]

In structural mechanics problems, K, C and M represent the \(n \times n\) stiffness, damping and mass matrices, respectively. a is the displacement, da ( = \(\dot{\mathbf{a}}\) ) is the velocity and d2a ( = \(\ddot{\mathbf{a}}\) ) is the acceleration.

The matrix f contains the time-dependent load vectors. If no external loads are active, use [] for f. The matrix f is organized as:

\[\begin{split}f = \begin{bmatrix} \text{time history of the load at } dof_1 \\ \text{time history of the load at } dof_2 \\ \vdots \\ \text{time history of the load at } dof_n \end{bmatrix}\end{split}\]

The dimension of f is:

(number of degrees-of-freedom) × (number of timesteps + 1)

The initial conditions are given by the vectors a0 and da0, containing initial displacements and initial velocities.

The matrix bc contains the time-dependent prescribed displacement. If no displacements are prescribed, use [] for bc. The matrix bc is organized as:

\[\begin{split}bc = \begin{bmatrix} dof_1 & \text{time history of the displacement} \\ dof_2 & \text{time history of the displacement} \\ \vdots & \vdots \\ dof_{m_2} & \text{time history of the displacement} \end{bmatrix}\end{split}\]

The dimension of bc is:

(number of dofs with prescribed displacement) × (number of timesteps + 2)

The time integration procedure is governed by the parameters given in the vector ip defined as:

\[ip = [dt, T, \alpha, \delta]\]

where dt specifies the time increment, T the total time, and alpha and delta are time integration constants for the Newmark family of methods.

Frequently used values:

\(\alpha\)	\(\delta\)	Method
\(\frac{1}{4}\)	\(\frac{1}{2}\)	Average acceleration (trapezoidal) rule
\(\frac{1}{6}\)	\(\frac{1}{2}\)	Linear acceleration
0	\(\frac{1}{2}\)	Central difference

The computed values of \(\mathbf{a}\), \(\dot{\mathbf{a}}\) and \(\ddot{\mathbf{a}}\) are stored in a, da and d2a, respectively. The first column contains the initial values and the following columns contain the values for each time step.

The dimension of a, da and d2a is:

(number of degrees-of-freedom) × (number of time steps + 1)

If the values are to be stored only for specific times, the parameter times specifies at which times the solution will be stored. The values are stored in a, da and d2a, one column for each requested time according to times. The dimension is then:

(number of degrees-of-freedom) × (number of requested times + 1)

If the history is to be stored in ahist, dahist and d2ahist for some degrees of freedom, the parameter dofs specifies for which degrees of freedom the history is to be stored. The computed time histories are stored in ahist, dahist and d2ahist, one row for each requested degree of freedom according to dofs. The dimension is:

(number of specified degrees of freedom) × (number of timesteps + 1)

In most cases only a few degrees-of-freedom are affected by the exterior load, and hence the matrix contains only few non-zero entries. In such cases it is possible to save space by defining f as sparse (a MATLAB built-in function).

sweep¶

Purpose:

Compute complex frequency response functions.

Syntax:

Y = sweep(K, C, M, p, w)

Description:

sweep computes the complex frequency response function for a system of the form:

\[[\mathbf{K} + i\omega\mathbf{C} - \omega^2 \mathbf{M} ]\mathbf{y}(\omega) = \mathbf{p}\]

Here, K, C, and M represent the m-by-m stiffness, damping, and mass matrices, respectively. The vector p defines the amplitude of the force. The frequency response function is computed for the values of \(\omega\) given by the vector w.

The complex frequency response function is stored in the matrix Y with dimension m-by-n, where n is equal to the number of circular frequencies defined in w.

Example:

The steady-state response can be computed by:

X = real(Y * exp(i * w * t));

and the amplitude by:

Z = abs(Y)