5.1.1. Newton’s equation of motion

Isaac Newton enunciated [NEW 87]:

The change of motion is proportional to the motive force impressed; and is made in the direction of the right fine in which that force is impressed.

Newton’s book is mainly geometrical as he did not have simple ways for the conceptualization of second derivatives at his disposition. We do not know if Newton abstracted acceleration as the second derivative of displacement with respect to time — this step implies the use of symbolic calculus and was performed by Jacob Hermann [HER 16] and, later, by Euler [EUL 50b]. Today, we write Newton’s second law as force = change in motion = change in [mass x velocity], i.e.

where F is the force applied to the system. Here, the force applied by the spring on the mass is − k(x − ℓ) = −k(u + x₀ − ℓ) and we have F = −k(u + x₀ − ℓ) + f. By assuming that x₀ = ℓ, this expression reduces to F = f − ku. When the mass is constant, we have

which leads to

i.e. force = mass x change in velocity = mass x acceleration.

5.1.2. Boundary value problem and initial condition problem

In the framework of Ordinary Differential Equations, it is usual to distinguish two kinds of differential equations: boundary value problems (BVP) and initial condition problems (ICP).

For a BVP, the values of the unknown are given on the extremities (boundaries) of the region of interest. For instance, if the analysis is performed on the time interval (0, T), we give the values of x or u at times t = 0 and t = T. For instance,

For ICP, we give the values at the left value of the interval (i.e. at t = 0):

As shown in the following, variational methods do not distinguish these two situations: both are treated by the same approach.

5.1.3. Generating a variational formulation

Variational formulations — also referred as variational equations — are obtained by performing three steps:

• 1) Multiplying the equation by a virtual motion y and integrating the result on the domain of interest (here, the interval (0, T)).
• 2) Using integration by parts (or Green’s Theorem, in the multidimensional case) in order to reduce the degree of derivation of x. Ideally, x and y must have the same degree in the final formulation, in order to get a symmetric expression.
• 3) Choosing the space V of the virtual motions. Ideally, the unknown x to be determined and the virtual motion y are both members of the same space — what may request some transformation on x — and the integrals correspond to scalar products on V. This step is not uniquely determined: different choices for V are possible and lead to different formulations. In this step, we must keep in mind that the differential equation is valid on the open domain Ω but not on its boundary ∂Ω: values on the boundary are to be considered as unknowns, except if they are given by a boundary condition.

For instance, let us consider the BVP described above.

In this case, Step 1 corresponds to write the equation as

and

Step 2 corresponds to the integration by parts

Step 3 is not uniquely determined. For instance, we may set

and we obtain a first variational formulation that reads as: find the unknowns P_T, P₀ and u V such that

As an alternative, we may set:

and

Then, a second variational formulation is given by: find w such that,

then determine .

Now, let us consider the ICP

In this case, Step 1 corresponds to

and

Step 2 corresponds yet to the integration by parts, but taking into account the initial conditions, we have :

Again, Step 3 is not uniquely determined. For instance, we may set

and we obtain a first variational formulation: find the unknowns P_T and u V such that

Alternatively, we may set

and

Then,

Alternatively, we may set

and

Then,

5.1.4. Generating an approximation of a variational equation

Analogously to algebraic equations, the solution is approximated as

with

The coefficients u_p , for p = 1, ... , k are unknowns to be determined. The unknowns may be ranged in a vector X including the coefficients and additional unknowns. X is solution of a linear system obtained by taking y(t) = φ_j(t) in the variational formulation and considering the boundary conditions that were not taken into account in the variational formulation.

5.1.5. Application to the first variational formulation of the BVP

For instance, let us consider the first variational formulation of the BVP: we consider

Then, the unknowns are X = (u₁, ... , u_k, P₀, P_T) and we have, for 1 ≤i ≤k,

This set of equations corresponds to a linear system A. X = B, where, for 1 ≤ i ≤ k:

For i = k + 1 and i = k + 2, we have:

5.1.6. Application to the second variational formulation of the BVP

Now, let us consider the alternative variational formulation of the BVP. In this case, we consider

and

Then, the unknowns are X = (w₁, ... , w_k) and, recalling that

we have

This set of equations corresponds to a linear system A. X = B, where, for 1 ≤ i ≤ k:

5.1.7. Application to the first variational formulation of the ICP

In this case, the unknowns to be determined are X = (u_{1, ... ,}u_nb, P_T) and we have, for 1 ≤ i ≤ k,

u(0) = u₀.

Thus, A. X = B, where, for 1 ≤ i ≤ k:

For i = k + 1:

5.1.8. Application to the second variational formulation of the ICP

In the second ICP formulation,

The unknowns to be determined are X = (u₁, ... , u_k, P_T) and we have for 1 ≤ i ≤ k,

Thus, A. X = B, where, for 1 ≤ i ≤ k:

5.2.1. From a variational equation to a differential equation

For an open domain Ω, the set (Ω)

plays a key role in the transformation. (Ω) is a set of virtual motions (or test functions) having the following useful properties:

• 1) Any φ (Ω) has infinitely many derivatives
• 2) For any φ (Ω): on the boundary of Ω, φ and all its derivatives are null.

The main result for the transformation of variational equations into differential equations is the following one:

THEOREM 5.1.— Let V ⊂ [L²(Ω)]ⁿ be a subspace such that [(Ω)]ⁿ ⊂V. Then

It admits the following extension

THEOREM 5.2.— Let V ⊂ [L²(Ω)]ⁿ be a subspace such that [(Ω)]ⁿ ⊂ V. Assume that Γ ⊂ ∂Ω, A [L²(Ω)]ⁿ, B [L²(Γ)]ⁿ, C_i ⁿ, x_i and

Let Then A = 0 on Ω − P, B = 0 on ∂Ω − Γ − P, C_i = 0 on P, (i = 1, ... , k).

In order to apply these theorems, we must eliminate all the derivatives of v from the variational equation. This is usually obtained by using either integration by parts:

or Green’s formula

EXAMPLE 5.1.— Let us consider {v H¹(Ω): v (0) = 0} and the variational equation

The first step is the elimination of the derivatives on v: we have

Thus,

and we have

Since (Ω) ⊂ V, the theorem shows that

Thus

EXAMPLE 5.2.— Let us consider and the variational equation

By using that

and Green’s formula, we have

Since v = 0 on Γ, we have

and

what yeilds

Since (Ω) ⊂ V, the theorem shows that

Thus,

5.2.2. From a differential equation to a variational equation

The converse transformation consists in the transformation of the differential equation into a Principle of Virtual Works. This is achieved by the multiplication by a virtual motion v and the subsequent use of integration by parts or Green’s formula. We are generally interested in obtaining a formulation where the virtual motion and the solution are members of the same space V. Initially, the space V is formed by all the functions defined on the domain of interest for which the quantities manipulated have a sense, but conditions satisfied by the solution and that are not used in order to generate the variational formulation have to be introduced in the definition of V.

EXAMPLE 5.3.— Let us consider the differential equation:

We start by the creation of a table containing:

• 1) the available equations;
• 2) the boundary or initial conditions;
• 3) the tools to be used: virtual works;
• 4) the steps to be performed:
  • i) integration by parts or use of Green’s formula;
  • ii) formal definition of the space.
• 5) The result of the steps performed.

For instance:

The aim is to use all the conditions. At the end of the transformation, any unused condition will be introduced in the definition of V.

As previously observed, the set of all the virtual motions is initially: V⁽⁰⁾ = {v. (0,1) → }. We start the transformation by using the differential equation itself:

In order to produce a virtual work, two steps are necessary: first, the equation is multiplied by an arbitrary virtual motion v:

Second, the resulting equality is integrated on Ω = (0,1):

Now, we use integration by parts:

Thus,we have

and

The condition u(0) = 0 cannot be used, since there is no u(0) in the variational formulation. However, the expression contains u′(1) and the condition u′(1) = 0 may be used:

The condition u(0) = 0 has not been used. Thus, it is introduced in the definition of the space:

Now, we have

The last step is the formal definition of the space: regarding to the expression, the first derivatives u′, v′ appear in a scalar product. Thus, a convenient choice is a space where the derivatives are square summable:

This is the final formulation: all the conditions were used.

EXAMPLE 5.4.— Let us consider the differential equation:

We have

Conditions and tools	Not used	Used
u = 0 on Γ	×
	×
f + Δu = 0 on Ω	×
Virtual work	×
Integration by parts or Green’s Formula	×
Formal definition of the space	×
Formulation:

We set V⁽⁰⁾ = {v: Ω → }.

Conditions and tools	Not used	Used
u = 0 on Γ	×
	×
f + Δu = 0 on Ω		✓
Virtual work	×
Integration by parts or Green’s Formula	×
Formal definition of the space	×
Formulation:

The equation is multiplied by an arbitrary virtual motion v:

The resulting equality is integrated on Ω:

Then, we use Green’s formula

Thus, we have

Conditions and tools	Not used	Used
u = 0 on Γ	×
	×
f + Δu = 0 on Ω		✓
Virtual work		✓
Integration by parts or Green’s Formula	×
Formal definition of the space	×
Formulation:

Since ∂u/∂n appears in the formulation, the corresponding condition may be used: we have

and

Conditions and tools	Not used	Used
u = 0 on Γ	×
		✓
f + Δu = 0 on Ω		✓
Virtual work		✓
Integration by parts or Green’s Formula
Formal definition of the space	×
Formulation:

The condition u = 0 on Γ has not been used. Thus, it is introduced in the definition of the space:

Now, we have

which corresponds to the following table:

Conditions and tools	Not used	Used
u = 0 on Γ		✓
		✓
f + Δu = 0 on Ω		✓
Virtual work		✓
Integration by parts or Green’s Formula		✓
Formal definition of the space	×
Formulation:

The formal definition of the space is obtained by considering the last formulation: gradients ∇u,∇v appear in a scalar product. Thus, analogously to the preceding situation, a convenient choice is a space where the derivatives are square summable:

Conditions and tools	Not used	Used
u = 0 on Γ		✓
		✓
f + Δu = 0 on Ω		✓
Virtual work		✓
Integration by parts or Green’s Formula		✓
Formal definition of the space		✓
Formulation:

This is the final formulation: all the conditions were used.

5.4.1. Linear equations

Let us assume that u → a(u, v) is linear. Using the approximation introduced in sections 4.4, 4.7 and 4.8:

and

where

In addition, we may approximate

In this case, we have

Thus,

These observations show that approximations of a linear differential equation using a total family or a Hilbert basis leads to the construction of:

• 1) a mass matrix M;
• 2) a stiffness matrix K;
• 3) a second member B, which may depends upon t;
• 4) an initial data U₀.

For instance, these matrices may be generated by finite element methods (FEM) or the continuous families presented.

The resulting equations may be integrated by intrinsic solvers from Matlab^® (ode45, ode23s, etc.). For instance, we may implement a subprogram evaluating M⁻¹(B − KU). Denoting by invM the inverse M⁻¹ of M and invMK the value of M⁻¹K, we may use the code

function v = second_member(t,U,invM,invMK, B)
v = invM*B(t) - invMK*U;
end

Then, we call one of these solvers:

T = ...; % give the time interval (0,T)
solver = @solvername; % solvername: ode45,ode23s,etc.
invM = M  eye(size(M)); % evaluate inv(M)
invMK = M  K; % evaluate inv(M)*K
U0= M  Y; % vector of initial conditions
B = @program_name; % name of the subprogram evaluating
F.
secm = @(t,X) second_member(t,U,invM,invMK, B);
sol = solver(secm, [0,T], U0);

Returning structure sol contains:

• 1) sol.x : vector of time moments. The line tc = sol.x; copies the time moments in vector tc — the solution has been evaluated at tc(1), tc(2) — tc(end).
• 2) sol.y : array of values of U. Each colum of sol.y contains the solution for one time moment: column j (sol. y (: , j) ) contains the solution U at time tc(j). Each line contains the solution for a given degree of freedom (dof) and all the times: line i (sol. y (i , :) ) contains the evolution of U_i for all the values tc.

5.4.2. Nonlinear equations

Let us assume that u → a(u, v) is nonlinear. In this case,

and we have

where

Let us assume that proj(U,basis) furnishes Pu, while a (pu, basis) furnishes a vector v such that v_i = a(PU, φ_i). In this case, the second member reads as

function v = second_member(t,U,invM,a, proj,basis,B)
pu = proj(U,basis);
aux = a(pu,basis) ;
v = invM*(B(t) -aux); % or M  (B(t) - aux) ;
end

Then, the differential equation may be solved by an intrinsic solver from Matlab^®.

5.4.3. Motion equations

Let us consider the equations describing the motion of an undamped system:

The model must be completed by initial conditions:

The main assumption on this system concerns the existence of normal modes, i.e. the existence of vectors φ_i, i = 1, … , n (n = dimension of the square matrices M, K and of vector U), such that the matrix Φ = [ϕ₁, …, ϕ_n] having as columns the vectors simultaneously diagonalizes the matrices M, K:

Where are diagonal matrices. The vectors are determined by setting

and solving

As an alternative, we may set

w is called a natural pulsation. Φ is the associated eigenvector, called normal mode.

We set

and the problem is rewritten as:

These equations correspond to evolution equations as treated in section 5.4.1. They may be integrated by intrinsic solvers from Matlab^® (ode45, ode23s, etc.). The subprogram evaluating the right-hand side of this system reads as

function v = secm_matlab(t,X,invM, invMK,invMC,F)
n = length(X)/2;
v = zeros(size(X));
v(1: n) = X(n+1:2*n);
v(n+1:2*n) = invM*F(t) - invMK*X(1: n) -
invMC*X(n+1:2*n);
end

Subprogram F must be implemented independently and returns the force vector at time t.

Assuming that the values of U₀, are defined in the Matlab vectors U0 and Up0, respectively, then Matlab^® solver is called as follows:

T = ...; % give the time interval (0,T)
solver = @solvername; % solvername: ode45,ode23s,etc.
invM = M  eye(size(M)); % evaluate inv(M)
invMK = M  K; % evaluate inv(M)*K
invMC = M  C; % evaluate inv(M)*C
X0 = [U0; Up0]; % vector of initial conditions
F = @program_name; % name of the subprogram evaluating
F.
secm = @(t,X) sec_matlab(t,X,invM, invMK,invMC, F);
sol = solver(secm, [0,T], X0);

5.4.3.1. Newmark integration of motion equations

One of the most popular numerical methods in dynamics is Newmark’s scheme:

Here, β, γ are strictly positive parameters to be defined by the user. For the method is unconditionally stable. Typical values are For the system under consideration, we have:

or, equivalently,

Let us introduce a time discretization involving nt subintervals of time:

We denote

Thus

The first equation reads as

and determines Uⁱ⁺¹. Then, the second equation determines . A simple Matlab implementation is the following:

nt = …; % number of time intervals
dt = T/nt; % time step
F = @(t) program_name(t, …); % subprogram evaluating F.
betta = …; % value of beta.
gama = … ; % value of gama.
A = MM + betta*dt*dt*KK;
t = 0;
Ui = U0;
Vi = Up0;
Fi = F(t);
Ai = M  (Fi - K*Ui)
for npt = 1: nptmax
t = t + dt;
Fi1 = F(t);
    aux = M*Ui + dt*M*Vi + 0.5*dt*dt*((1-2*betta)*(Fi -
K*Ui) + 2*betta*Fi1);
    Ui1 = A aux;
    Ai1 = M  (Fi1 - K*Ui1);
Vi1 = Vi + dt*((1-gama)*Ai + gama*Ai1);
Ui = Ui1;
    Vi = Vi1;
    Ai = Ai1;
Fi = Fi1 ;
    end;

Newmark’s scheme also reads as the construction of a sequence of vectors given by

or

Since

we have

so that Newmark’s scheme reads as

with

The first equation furnishes a linear system for U^k+1. Once U^k+1 is determined, the second line furnishes . An alternative code is

nt = …; % number of time intervals
dt = T/nt; % time step
F = @program_name; % name of the subprogram evaluating
F.
betta = …; % value of beta.
gama = … ; % value of gama.
A1_11 =M + betta*dt*dt*K;
A1_21 = gama*dt*K;
A1_22 = M;
A2_11 = M - (0.5 - betta)*dt*dt*K;
A2_12 = dt*M;
A2_21 = -(1-gama)*dt*K;
A2_22 = M;
A3_11 = (0.5 - betta)*dt*dt;
A3_12 = betta*dt*dt ;
A3_21 = (1-gama)*dt;
A3_22 = (1-gama)*dt;
t = 0;
Uk = U0;
Upk = Up0;
Fk = F(t);
for npt = 1: nt
    t = t + dt;
Fk1 = F(t);
    B = A2_11*Uk + dt*M*Upk + A3_11*Fk + A3_12*Fk1;
    Uk1 = A1_11  B;
    B = - gama*dt*K*Uk1-(1-
gama)*dt*K*Uk+M*Upk+A3_21*Fk+A3_22*Fk1;
    Upk1 = M  B;
    Uk = Uk1;
    Upk = Upk1;
    Fk = Fk1;
end;

5.4.3.2. Fourier integration of motion equations

If

then

with

Here, ϕ_k are the eigenvectors and w_k,nat are the natural pulsations associated to K. We have

The coefficients , UA, UB are determined by

Let us introduce the matrix Φ = [ϕ₁, … ,ϕ_n] having as columns the eigenvectors, Ω_nat a diagonal matrix containing the natural frequencies in the diagonal, SIN_nat and COS_nat, diagonal matrices containing the terms sin(w_nat,kt) , cos(w_nat,kt), respectively:

Then:

and

Analogously, let us consider Ω a diagonal matrix containing the frequencies w_i in the diagonal, SIN and COS, diagonal matrices containing the terms sin(w_it) , cos(w_it), respectively:

Then

so that

Thus

which determines the coefficients A, B. A simple implementation of this method is the following one:

function [Ubar, U_A,U_B,a,b] =
solve_multichromatic(omega, M,K, Phi, om_nat, Fbar, F_A,
F_B, U0, V0)
U_A = cell(size(omega));
U_B = cell(size(omega));
Ubar = K Fbar;
for i = 1: length(omega)
om = omega(i);
    fA = F_A{i};
    fB = F_B{i};
    matsys = K - om^2*M;
uA = matsys  fA;
    uB = matsys  fB;
U_A{i} = uA;
    U_B{i} = uB;
end;
X_B = Ubar;
for i = 1:length(omega)
X_B = X_B + U_B{i};
end;
b = Phi  (U0 - X_B);
X_A = zeros(size(U_A{1}));
for i = 1:length(omega)
    X_A = X_A + omega(i)*U_A{i};
end;
gama = Phi  (V0 - X_A);
a = gama./om_nat;
return;
end

5.4.4. Motion of a bar

Let us consider a bar of length ℓ > 0, submitted to a distributed load f , clamped at origin and force f_ℓ extremity ℓ. The motion’s equations are

where:

• 1) E = Young’s modulus of the material;
• 2) S = surface of the right section;
• 3) ρ = volumetric density of the material.

FEM discretization starts by the generation of a set of vertices and elements, defined by a partition

Element number i is (x_i, x_i+1). When a uniform discretization is used, we have

In this case, we may generate the discretization and the elements by using the following matlab^® code

ele = …; % give the value of ℓ
n = …; % give the number of elements
x = linspace(0, ele, n+1);
element = zeros(n,2);
for i = 1: n
    element(i,:) = [i i+1];
end;

The complete mass and stiffness cM, cK are generated by using integrals involving the shape functions on the elements, i.e. the elementary matrices m_e for mass and k_e for stiffness. For instance, let us consider the classical bar element involving 2 degrees of freedom by element. In this case,

For instance, we may use the Matlab code

   elastmod = …; % give the value of Young's modulus
   surf = …; % give the value of the surface of the
right section
   rho = … ; % give the value of the volumetric density
of the material
   m_e = rho*surf*h*[2 1 ; 1 2]/6;
   k_e = elastmod*surf*[1 -1; -1 1]/h;

Matrices cM, cK are obtained by assembling the elementary matrices on the elements:

cM = zeros(n+1,n+1);
cK = zeros(n+1,n+1);
for i = 1:n
   ind = element(i,:);
   cM(ind, ind) = cM(ind, ind) + m_e;
   cK(ind, ind) = cK(ind, ind) + k_e;
end;

Boundary condition u(0, t) = 0 is taken into account by the elimination of the first line and column of cM, cK. For instance, we may use the code:

M = cM(2:n+1,2:n+1);
K = cK(2:n+1,2:n+1);

External forces are given by a time-dependent field f(x, t). Force vector is obtained by assembling the elementary forces: on element i, corresponding to nodes [i i+1] and degrees of freedom [i i+1], the forces are (f_i, f_i+1), f_i = f(x_i, t). Assembling is made by using the vector of nodal forces f = (f₁, ... f_n+1)^t. Let us assume that the values of f, f_ℓ are furnished by subprograms f, f_ele, respectively. Then, an example of assembling is given by the subprogram

   function cF = assembly_force(f,f_ele,x,n,t, elastmod,
surf)
   cF = zeros(n+1,1);
   for i = 1: n
   dof = [i, i+1];
       f_el = [f(x(dof(1)),t), f(x(dof(2)),t)];
       cF(dof) = cF(dof) + mean(f_el)*[1; 1];
end;
cF(n+1) = CF(n+1) + f_ele(t)/(elastmod*surf);
return;
end

In this situation, assembling may be simplified:

At time t, this vector may be obtained by using the subprogram

   function cF = assembly_force(f,f_ele, x,n,t, cM,
elastmod, surf)
   f = zeros(n+1,1);
   for i = 1: n+1
   f(i) = f(x(i),t);
   end;
   cF = cM*f;
   cF(n+1) = cF(n+1) + f_ele(t)/(elastmod*surf); );
   return;
   end

The elimination of the first degree of freedom is performed as follows:

   function F = force(t,f,f_ele, x,n, cM, elastmod, surf)
   cF = assembly_force(f,f_ele, x,n,t, cM, elastmod,
surf);
   F = cF(2:end);
   return;
   end

Initial conditions are discretized as follows: assume that the values of u₀(x), u₀(x) are furnished by subprograms u_0, up_0, respectively. Then, we may use the code

cU0 = zeros(n+1,1);
cUp0 = zeros(n+1,1);
for i = 1: n+1
    cU0(i) = u_0(x(i));
cUp0(i) = up_0(x(i));
end;
U0 = cU0(2:n+1);
Up0 = cUp0(2:n+1);

EXAMPLE 5.7.— A uniform monochromatic force

Let us consider the bar submitted to the oscillating forces

For

omega = 1;
pA = 1;
pB = 2;
f = @(x,t) pA*sin(omega*t) + pB*cos(omega*t);
f_ele = @(t) 0;

we have the results below for the motion of the last node.

EXAMPLE 5.8.— A uniform multichromatic force

Let us consider the bar submitted to the oscillating forces

For

omega = [ 1 2 3 4];
pA = [ 4 3 2 1];
pB = [ 1 2 3 4];
pbar = 0;
f = @(x,t) somme_trigo(x,t, pbar, pA, pB, omega);
f_ele = @(t) 0;

we have the results below for the motion of the last node.

5.4.5. Motion of a beam under flexion

Let us consider a beam of length ℓ > 0, submitted to pure flexion under a distributed load f, clamped at origin and force f_ℓ, extremity ℓ. The motion’s equations are

where:

• 1) E = Young’s modulus of the material;
• 2) S = surface of the right section;
• 3) I = inertia moment of the cross section;
• 4) ρ = volumetric density of the material.

FEM discretization is analogous to the preceding one. For instance, let us consider the classical beam element involving 4 degrees of freedom by element. In this case,

For instance, we may use the Matlab code

   elastmod = ...; % give the value of Young's modulus
   surf = ...; % give the value of the surface of the
right section
   rho = ... ; % give the value of the volumetric density
of the material
   inertia = ... ; % give the value of the inertia of the
right section
   m_el = rho*surf*h*[
   156 22*h 54 -13*h;
              22*h 4*h^2 13*h -3*h^2;
              54 13*h 156 -22*h;
              -13*h -3*h^2 -22*h 4*h^2
          ]/420;
   k_el = elastmod*inertia*[
          12 6*h -12 6*h;
          6*h 4*h^2 -6*h 2*h^2;
          -12 -6*h 12 -6*h;
          6*h 2*h^2 -6*h 4*h^2
   ]/h^3;

Analogously to the preceding situation, matrices cM, cK are obtained by assembling the elementary matrices on the elements:

ndof = 2*(n+1);
cM = zeros(ndof,ndof);
cK = zeros(ndof,ndof);
for i = 1: n
ind = (2*i-1:2*i+2);
    cM(ind,ind) = cM(ind,ind) + m_el;
    cK(ind,ind) = cK(ind,ind) + k_el;
end;

Boundary conditions u(0, t) = u′(0, t) = 0 are taken into account by the elimination of the two first lines and columns of cM, cK. For instance, we may use the code:

M = cM(3:ndof,3:ndof);
K = cK(3:ndof,3:ndof);

External forces are given by a time-dependent field f(x, t). The force vector is obtained by assembling the elementary forces: on element i, corresponding to nodes [i i+1] and degrees of freedom [2*i−1 2*i 2*i+1 2*i+2], the forces are (f_i, f_i+1), f_i = f(x_i, t). An example of assembling is given by the subprogram

   function cF = assembly_force(f,f_ele,x,n,t, elastmod,
inertia)
   cF = zeros(2*n+2,1);
   for i = 1: n
       dof = (2*i-1:2*i+2);
       f_el = [f(x(i),t); f(x(i+1),t)];
       cF(dof) = cF(dof) + mean(f_el)*h*[1; h/6; 1; -
h/6]/2;
   end;
   cF(2*n+1) = CF(2*n+1) + f_ele(t)/(elastmod*inertia);
   return;
   end

The elimination of the two first degrees of freedom is performed as follows:

   function F = force(t,f,f_ele, x,n, cM, elastmod,
inertia)
   cF = assembly_force(f,f_ele, x,n,t, cM, elastmod,
inertia);
   F = cF(3:end);
   return;
   end        

Initial conditions are discretized analogously to the case of a bar.

EXAMPLE 5.9.— A uniform monochromatic force

Let us consider the beam submitted to the oscillating forces

For

omega = 1;
pA = 1;
pB = 2;
f = @(x,t) pA*sin(omega*t) + pB*cos(omega*t);
f_ele = @(t) 0;

we have the results below for the motion of the last node.

EXAMPLE 5.10.— A uniform multichromatic force

Let us consider the beam submitted to the oscillating forces

For

omega = [ 1 2 3 4];
pA = [ 4 3 2 1];
pB = [ 1 2 3 4];
pbar = 0;
f = @(x,t) somme_trigo(x,t, pbar, pA, pB, omega);
f_ele = @(t) 0;

we have the results below for the motion of the last node.