Christmas Calculation - the beauty of Euler's Formula

How a structure responds to an earthquake

Following my article on earthquakes, a few of my technically savvy readers, and structural engineers in particular, have asked for more details on the relationship between structures and control systems. For those readers without a background in mathematics, I apologize because this article might be a bit too complex to follow.

However, for those of us who have taken a course in differential equations, we might recall Euler’s formula, which the physicist Richard Feynman referred to as 'the most remarkable formula in mathematics”.

\(\begin{align} e^{it} &= \cos(t) + i \sin(t) \\ \end{align}\)

Let’s construct a parametric function using this and plot it in the complex plane. Although it is not always stated, the above expression can be reduced to a function of four dimensions: Z, X, Y, and t, with t being the chosen independent parameter. For convenience, it is often thought as an angle, but it’s worth noting that all angles are essentially just fractions - that can be expressed as rational numbers.

\(\begin{align} x &= \cos(t) \\ y &= i\sin(t) \\ z &= x + y \end{align}\)

If I were to plot three of the above variables, noting that y is complex, I would obtain the following graph, which shows that Euler’s formula is actually the equation of a spring.

The distance of the spring from the origin and the tangent angle are determined using the following expression.

\(|z| = \sqrt{x^2 + y^2}, \quad \theta = \tan^{-1}\left( \frac{y}{x} \right) \)

Because of the nature of Euler's formula, the spring's projection onto the complex plane forms a circle. In the other planes, it forms classical cosine and sine waves.

Up until this point, there is only mathematics, but now recall that I mentioned earthquakes are modeled using an extension of Newton’s second Law. The basic equation of motion for a harmonic oscillating spring is given by the following expression.

\(ma +cv + kx = 0 \)

The first is D'Alembert's principle of inertial and the last is Hooke’s law for an elastic spring. The second term represents friction, often referred to as damping. If it were neglected, as is the case in the above graph, the spring would vibrate indefinitely at its natural frequencies and the undamped homogeneous solution would be obtained.

If I were to express velocity and acceleration it in differential terms and add the earthquake force p(t) as a boundary condition, I would obtain the following expression:

The traditional approach by structural engineers is to solve the above equation using Duhamel’s Integral. However, it can just as easily be solved using Laplace transformation, resulting in the following expression if p(t) = ma_g, with a_g the acceleration of the ground.

\(G(s) = \frac{X(s)}{F(s)} = \frac{1}{ms^2 + cs + k}\)

if s = i*omega, the solution would reduce to the Fourier transform.

By now, electrical engineers should understand what is happening, but I have found that many structural engineers tend to become lost. This is because we often think in terms of bending moments and shear forces, without realizing that they are also solutions to differential equations.

The dynamic behavior of a structure can be explained as a set of coupled springs, which can be regarded as a system that “filters” the noise. The behavior of the system is described by a second-order differential equation. The distance between any two points is considered the transfer function that describes the filter.

These equations are comparable to those that electrical engineers refer to as control systems.

A structure is a control system, and in fact if you were to solve your moment-curvature equations, such as the Euler-Bernoulli Beam, using Laplace or Fourier Transform, the solution would would converge to the classical solutions in our textbooks. The results would be the exact same as those that we obtain from Bending Moment and Shear Force Diagrams or those that we obtain from deflection curves.

Think of it this way, our structures "TRANSFER” loads from one point to another, hence they can be described by transfer functions. Therefore they are also control systems.

For those who want to play around with it, I added a simple python script that is below this letter.

"""
By Hügo KRÜGER
This code visualizes a 3D projection for a function Z, where:
t - Input parameter (time)
X - Real part of Z
Y - Imaginary part of Z
The 3D plot shows projections along T (time), X (real), and Y (imaginary) axes.
Magnitude is determined as abs(Z), and phase is calculated as tan(Y/X).
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generate time range and the complex function Z
t = np.linspace(-3, 3, 1000)
Z = np.exp(t * np.pi * 1j)

# Extract real and imaginary parts, and compute the phase in degrees
real_part = np.real(Z)
imag_part = np.imag(Z)
phase_deg = np.degrees(np.angle(Z))

# Set axis limits for symmetry based on real and imaginary parts
axis_limit = max(abs(real_part).max(), abs(imag_part).max())

# Create grids for projections on planes
mesh_size = 50
xy_mesh = np.linspace(-axis_limit, axis_limit, mesh_size)
time_mesh = np.linspace(min(t), max(t), mesh_size)

# 3D Plot setup
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# Set axis labels, limits, and title
ax.set_xlim([-axis_limit, axis_limit])
ax.set_ylim([-axis_limit, axis_limit])
ax.set_zlim([min(t), max(t)])
ax.set_xlabel('Real Part')
ax.set_ylabel('Imaginary Part')
ax.set_zlabel('Time (t)')
ax.set_title('3D Plot of Z and Projections')

# Plot projection planes
# XY-plane (Real vs Imaginary)
X, Y = np.meshgrid(xy_mesh, xy_mesh)
ax.plot_surface(X, Y, np.zeros_like(X), color='gray', alpha=0.3)

# Time-X plane (Time vs Real Part)
X2, Y2 = np.meshgrid(xy_mesh, time_mesh)
ax.plot_surface(X2, np.zeros_like(X2), Y2, color='blue', alpha=0.3)

# Time-Y plane (Time vs Imaginary Part)
X3, Y3 = np.meshgrid(time_mesh, xy_mesh)
ax.plot_surface(np.zeros_like(Y3), Y3, X3, color='green', alpha=0.3)

# Scatter plot of Z colored by phase
sc = ax.scatter(real_part, imag_part, t, c=phase_deg, cmap='hsv', marker='o')

# Add projections of Z onto planes
ax.plot(real_part, imag_part, np.zeros_like(t), 'r-', alpha=0.7)  # Real-Imaginary
ax.plot(real_part, np.zeros_like(t), t, 'black', alpha=0.7)  # Real-Time
ax.plot(np.zeros_like(t), imag_part, t, 'b-', alpha=0.7)  # Imaginary-Time

# Add a colorbar for the phase
plt.colorbar(sc, ax=ax, label='Phase (Degrees)').set_ticks([-180, -90, 0, 90, 180])

# Show the 3D plot
plt.show()

# 2D Projections
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# 2D Real vs Imaginary plot
axes[0].plot(real_part, imag_part, 'r-', alpha=0.7)
axes[0].set_xlim([-axis_limit, axis_limit])
axes[0].set_ylim([-axis_limit, axis_limit])
axes[0].set_xlabel('Real Part')
axes[0].set_ylabel('Imaginary Part')
axes[0].set_title('Real-Imaginary Plane')
axes[0].grid(True)
axes[0].axis('equal')

# 2D Real vs Time plot
axes[1].plot(t, real_part, 'black', alpha=0.7)
axes[1].set_xlim([min(t), max(t)])
axes[1].set_ylim([-axis_limit, axis_limit])
axes[1].set_xlabel('Time (t)')
axes[1].set_ylabel('Real Part')
axes[1].set_title('Real-Time Plane')
axes[1].grid(True)

# 2D Imaginary vs Time plot
axes[2].plot(t, imag_part, 'b-', alpha=0.7)
axes[2].set_xlim([min(t), max(t)])
axes[2].set_ylim([-axis_limit, axis_limit])
axes[2].set_xlabel('Time (t)')
axes[2].set_ylabel('Imaginary Part')
axes[2].set_title('Imaginary-Time Plane')
axes[2].grid(True)

# Adjust layout and show the 2D plots
plt.tight_layout()
plt.show()