
From what we've seen so far, solving the Duhamel integral analytically is certainly possible, but as we've already mentioned, real-life engineering is messy and often we can't neatly describe our applied loading in a mathematical function. This is where the numerical solution strategy comes in. In this section we'll build a function that numerically solves the Duhamel integral for any arbitrary loading we pass into the function. We'll start by restating the Duhamel integral for a damped system,

$$
u(t) = \frac{1}{m\omega_d}\int_0^t F(\tau) e^{-\xi\omega_n(t-\tau)}\sin \omega_d (t-\tau)\mathrm{d}\tau
\tag{1}
$$

We can separate out the time variables with the help of a trigonometric identity,

$$
u(t) = \frac{1}{m\omega_d} e^{-\xi\omega_n t}\int_0^t F(\tau) e^{\xi\omega_n\tau} \underbrace{\left[ \sin (\omega_d t) \cos (\omega_d \tau) - \cos (\omega_d t) \sin (\omega_d \tau)  \right]}_{\text{expanded portion}} \mathrm{d}\tau
\tag{2}
$$

Since we're integrating with respect to $\tau$ we can restructure our expression to isolate the terms requiring numerical integration,

$$
u(t) = A e^{-\xi \omega_n t} \sin (\omega_d t) - Be^{-\xi\omega_n t} \cos(\omega_d t)
\tag{3}
$$

where,

$$
A = \frac{1}{m\omega_d} \int_0^t e^{\xi \omega_n \tau} F(\tau)\cos(\omega_d\tau)\mathrm{d}\tau
\tag{4}
$$

$$
B = \frac{1}{m\omega_d} \int_0^t e^{\xi \omega_n \tau} F(\tau)\sin(\omega_d\tau)\mathrm{d}\tau
\tag{5}
$$

$A$ and $B$ now
represent two integrals that we can numerically evaluate. We'll use the trapezoidal
rule here; we've covered this in some detail in our [beam deflection calculator project](/project/building-a-beam-deflection-calculator) so we'll only briefly summarise the key points here.

If we consider the integrand (of $A$ or $B$) as a simple function, $f(\tau)$, we could plot this function on a $\tau-f(\tau)$ plane. We know that the integration is simply the area bounded by the function and the horizontal axis. Using the trapezoidal rule we can approximate this as,

$$
\text{Area} = \int_0^{t_n} f(\tau) \mathrm{d}\tau \approx \frac{\Delta_t}{2}\left[f(t_0) + 2f(t_1)+\ldots + 2f(t_{n-1}) + f(t_n)\right]
\tag{6}
$$

Now all that remains is to write a function to cycle through the time history computing the numerical integrations in equations 4 and 5 and then compute the response using equation 3.

## Implementing the numerical solution

We'll package all of our code into a function that accepts a time and force vector. This way we can easily reuse the function from this file and deploy it elsewhere down the road. Let's start by defining SDoF system parameters; mass $m$, damping ratio $\xi$ and natural frequency/angular natural frequency $f_n/\omega_n$. We can also define the damped natural frequency, $\omega_d = \omega_n\sqrt{1-\xi^2}$.

```python
m = 1000 #(kg) Mass
xi = 0.05 # Damping ratio
f = 1.5 #(Hz) Natural frequency
wn = 2*math.pi*f #(rads/s) Angular natural frequency
wd = wn*math.sqrt(1-xi**2) #(rads/s) Damped angular natural frequency
```

For the purposes of demonstration here we'll define a harmonic forcing function. This will make it easier to validate our Duhamel integral function against a closed-form solution. First we can define a time vector in the usual way.

```python
tMax = 20 #(s) Max time
delT = 0.01 #(s) Time-step
time = np.arange(0, tMax+delT, delT) #Time vector

f_Force = 1 #(Hz) Forcing frequency
wf = 2*math.pi*f_Force #(rads/s) Angular forcing frequency
P=100 #(N) Force amplitude
force = P*np.sin(wf*time) #Force vector
```

The function we write next will accept the time and force vectors as arguments. We'll rely on the fact that all of the system parameters are defined globally so we won't bother explicitly passing them into the function. We basically need to cycle through the time record and compute the area of the $i^{th}$ trapezoid and add this to a cumulative sum of areas as we move along the time axis. In this way we can calculate appropriate values for $A$ and $B$ (equations 4 and 5) at each time-step. Then, we can calculate the actual response for each time-step using equation 3. If I've lost you with the trapezoidal method of numerical integration, go back to the beam deflection tutorial I referenced earlier. I walk through this technique from scratch there.

```python
def Duhamel(T, F):
#Initialise a container of zeros to hold the displacement values
U = np.zeros(len(T))

#Initialise values for the cumulative sum used to calculate A and B at each time-step
ACum_i=0
BCum_i=0

#Cycle through the time vector and evaluate the response at each time point
for i, t in enumerate(T):
    #Only start calculating on the second iteration (need two values for trapezoidal area calculation)
    if i>0:
        #Calculate A[i] - equation 4
        y_i = math.e**(xi*wn*T[i]) * F[i] * math.cos(wd*T[i]) #Value at current time-step
        y_im1 = math.e**(xi*wn*T[i-1]) * F[i-1] * math.cos(wd*T[i-1]) #Value at previous time-step
        Area_i = 0.5*delT*(y_i+y_im1) #Area of trapezoid
        ACum_i += Area_i #Cumulative area from t=0 to current time
        A_i = (1/(m*wd))*ACum_i #Value of A for the current time-step

        #Calculate B[i] - equation 5 (same notes as for A above)
        y_i = math.e**(xi*wn*T[i]) * F[i] * math.sin(wd*T[i])
        y_im1 = math.e**(xi*wn*T[i-1]) * F[i-1] * math.sin(wd*T[i-1])
        Area_i = 0.5*delT*(y_i+y_im1)
        BCum_i += Area_i
        B_i = (1/(m*wd))*BCum_i

        #Calculate the response - equation 3
        U[i] = A_i*math.e**(-xi*wn*T[i])*math.sin(wd*T[i]) - B_i * math.e**(-xi*wn*T[i])*math.cos(wd*T[i])

return U
```

Now that we have the function defined we just need to call it and assign the output to a variable for plotting.

```python
response = Duhamel(time, force) #Response calculated using the Duhamel integral function

#Initialise a fig to plot onto
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])
axes.plot(time,response)

#Housekeeping
axes.set_xlabel('time (s)')
axes.set_ylabel('Disp (m)')
axes.set_title('SDoF system response to harmonic loading')
axes.set_xlim([0,tMax])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/numerical-solution/img1.jpg"
	alt="Response to harmonic loading. | EngineeringSkills.com"
	width="1000"
	height="369"
	caption="Fig 1. Response to harmonic loading."
/>

We can see above that the response contains a transient component that dissipates due to damping, leaving a steady-state response which is exactly what we would expect from the harmonic excitation of a SDoF system.


Please log in to continue