结构动力学仿真-主题043-结构动力学数值方法

kkchenkx

548人浏览 · 2026-04-06 21:05:29

kkchenkx · 2026-04-06 21:05:29 发布

主题043：结构动力学数值方法

1. 引言

1.1 数值方法的重要性

在实际工程问题中，结构动力学方程往往难以获得解析解。复杂的边界条件、非线性材料特性、不规则几何形状等因素使得解析方法无能为力。数值方法成为求解结构动力学问题的主要手段。

数值方法的核心思想是将连续的微分方程离散化为代数方程组，通过计算机求解获得近似解。随着计算机技术的发展，数值方法在结构动力学领域的应用越来越广泛，能够处理大规模、高复杂度的工程问题。
在这里插入图片描述

1.2 常用数值方法概述

结构动力学数值方法主要分为以下几类：

时间积分方法：

中心差分法（显式）
Newmark-β法（隐式）
Wilson-θ法
Houbolt法

空间离散方法：

有限差分法（FDM）
有限元法（FEM）
边界元法（BEM）
谱方法

模态分析方法：

子空间迭代法
Lanczos方法
Ritz向量法

模型降阶方法：

模态截断法
静力凝聚法
动态子结构法（Craig-Bampton）

1.3 数值方法的评价标准

评价一个数值方法的优劣通常考虑以下指标：

精度：数值解与精确解的接近程度
稳定性：长时间积分是否保持有界
收敛性：网格加密或步长减小时是否收敛到精确解
计算效率：计算时间和内存消耗
适用性：能否处理各种复杂问题

2. 有限差分法

2.1 有限差分法基本原理

有限差分法（Finite Difference Method, FDM）是最古老的数值方法之一，其基本思想是用差商代替导数，将微分方程转化为差分方程。

2.1.1 差分格式

对于函数 $u (t)$ ，常用的差分格式包括：

向前差分（一阶精度）：

$dudt∣tn≈un+1−unΔt\frac{du}{dt}\bigg|_{t_n} \approx \frac{u_{n+1} - u_n}{\Delta t}$

向后差分（一阶精度）：

$dudt∣tn≈un−un−1Δt\frac{du}{dt}\bigg|_{t_n} \approx \frac{u_n - u_{n-1}}{\Delta t}$

中心差分（二阶精度）：

$dudt∣tn≈un+1−un−12Δt\frac{du}{dt}\bigg|_{t_n} \approx \frac{u_{n+1} - u_{n-1}}{2\Delta t}$

二阶中心差分：

$d2udt2∣tn≈un+1−2un+un−1Δt2\frac{d^2u}{dt^2}\bigg|_{t_n} \approx \frac{u_{n+1} - 2u_n + u_{n-1}}{\Delta t^2}$

2.1.2 截断误差分析

利用Taylor展开可以分析差分格式的精度：

$\Delta t) = u(t) + \Delta t \cdot u'(t) + \frac{\Delta t^2}{2} u''(t) + \frac{\Delta t^3}{6} u'''(t) + O(\Delta t^4)$

$\Delta t) = u(t) - \Delta t \cdot u'(t) + \frac{\Delta t^2}{2} u''(t) - \frac{\Delta t^3}{6} u'''(t) + O(\Delta t^4)$

将两式相减：

$\Delta t) - u(t - \Delta t) = 2\Delta t \cdot u'(t) + \frac{\Delta t^3}{3} u'''(t) + O(\Delta t^5)$

因此：

$u(t+Δt)−u(t−Δt)2Δt=u′(t)+Δt26u′′′(t)+O(Δt4)\frac{u(t + \Delta t) - u(t - \Delta t)}{2\Delta t} = u'(t) + \frac{\Delta t^2}{6} u'''(t) + O(\Delta t^4)$

截断误差为 $O(Δt2)O(\Delta t^2)$ ，即二阶精度。

2.2 中心差分法求解运动方程

2.2.1 算法推导

考虑单自由度系统的运动方程：

$mu¨+cu˙+ku=F(t)m\ddot{u} + c\dot{u} + ku = F(t)$

用中心差分近似速度和加速度：

$u˙n=un+1−un−12Δt\dot{u}_n = \frac{u_{n+1} - u_{n-1}}{2\Delta t}$

$u¨n=un+1−2un+un−1Δt2\ddot{u}_n = \frac{u_{n+1} - 2u_n + u_{n-1}}{\Delta t^2}$

代入运动方程：

$mun+1−2un+un−1Δt2+cun+1−un−12Δt+kun=Fnm\frac{u_{n+1} - 2u_n + u_{n-1}}{\Delta t^2} + c\frac{u_{n+1} - u_{n-1}}{2\Delta t} + ku_n = F_n$

整理得到递推公式：

$(mΔt2+c2Δt)un+1=Fn−(k−2mΔt2)un−(mΔt2−c2Δt)un−1\left(\frac{m}{\Delta t^2} + \frac{c}{2\Delta t}\right)u_{n+1} = F_n - \left(k - \frac{2m}{\Delta t^2}\right)u_n - \left(\frac{m}{\Delta t^2} - \frac{c}{2\Delta t}\right)u_{n-1}$

令：

$a0=mΔt2+c2Δta_0 = \frac{m}{\Delta t^2} + \frac{c}{2\Delta t}$

$a1=k−2mΔt2a_1 = k - \frac{2m}{\Delta t^2}$

$a2=mΔt2−c2Δta_2 = \frac{m}{\Delta t^2} - \frac{c}{2\Delta t}$

则：

$un+1=Fn−a1un−a2un−1a0u_{n+1} = \frac{F_n - a_1 u_n - a_2 u_{n-1}}{a_0}$

2.2.2 启动计算

中心差分法是三步格式，需要 $u_0$ 和 $u_{-1}$ 来启动。 $u_0$ 由初始条件给出， $u_{-1}$ 可通过初始速度近似：

$u˙0=u1−u−12Δt\dot{u}_0 = \frac{u_1 - u_{-1}}{2\Delta t}$

因此：

$u−1=u1−2Δt⋅u˙0u_{-1} = u_1 - 2\Delta t \cdot \dot{u}_0$

或者利用Taylor展开：

$u−1=u0−Δt⋅u˙0+Δt22⋅u¨0u_{-1} = u_0 - \Delta t \cdot \dot{u}_0 + \frac{\Delta t^2}{2} \cdot \ddot{u}_0$

其中 $u¨0\ddot{u}_0$ 可由初始时刻的运动方程求得：

$u¨0=F0−cu˙0−ku0m\ddot{u}_0 = \frac{F_0 - c\dot{u}_0 - ku_0}{m}$

2.2.3 稳定性条件

中心差分法是条件稳定的，稳定性条件为：

$Δt≤Tnπ=2ωn\Delta t \leq \frac{T_n}{\pi} = \frac{2}{\omega_n}$

其中 $T_n$ 是系统的最小周期， $ωn\omega_n$ 是最高固有频率。

对于多自由度系统，稳定性条件为：

$Δt≤2ωmax\Delta t \leq \frac{2}{\omega_{max}}$

其中 $ωmax\omega_{max}$ 是系统的最高固有频率。

2.3 有限差分法的优缺点

优点：

算法简单，易于实现
计算效率高（显式格式，无需解方程组）
适用于大规模问题

缺点：

条件稳定，时间步长受限
对复杂几何和边界条件适应性差
精度相对较低

3. Newmark-β法

3.1 Newmark-β法基本原理

Newmark-β法是一类重要的隐式时间积分方法，由N.M. Newmark于1959年提出。该方法通过引入两个参数 $β\beta$ 和 $γ\gamma$ 来控制算法的精度和稳定性。

3.1.1 基本假设

Newmark-β法假设在时间步长 $t_n, t_{n+1}]$ 内，加速度和速度的变化规律为：

$u˙n+1=u˙n+(1−γ)Δt⋅u¨n+γΔt⋅u¨n+1\dot{u}_{n+1} = \dot{u}_n + (1-\gamma)\Delta t \cdot \ddot{u}_n + \gamma\Delta t \cdot \ddot{u}_{n+1}$

$un+1=un+Δt⋅u˙n+(12−β)Δt2⋅u¨n+βΔt2⋅u¨n+1u_{n+1} = u_n + \Delta t \cdot \dot{u}_n + \left(\frac{1}{2}-\beta\right)\Delta t^2 \cdot \ddot{u}_n + \beta\Delta t^2 \cdot \ddot{u}_{n+1}$

其中：

$γ\gamma$ 控制数值阻尼，通常取0.5（无人工阻尼）
$β\beta$ 控制精度，通常取0.25（平均加速度法）

3.1.2 常用参数组合

平均加速度法（无条件稳定，二阶精度）：

$γ=0.5\gamma = 0.5$
$β=0.25\beta = 0.25$

线性加速度法（条件稳定，二阶精度）：

$γ=0.5\gamma = 0.5$
$β=1/6\beta = 1/6$
稳定性条件： $Δt≤0.551Tn\Delta t \leq 0.551T_n$

Fox-Goodwin法：

$γ=0.5\gamma = 0.5$
$β=1/12\beta = 1/12$

3.2 Newmark-β法算法实现

3.2.1 算法步骤

对于线性系统 $Mu¨+Cu˙+Ku=F(t)M\ddot{u} + C\dot{u} + Ku = F(t)$ ，Newmark-β法的计算步骤如下：

1. 初始化

给定初始条件 $u_0$ 和 $u˙0\dot{u}_0$
计算初始加速度： $u¨0=M−1(F0−Cu˙0−Ku0)\ddot{u}_0 = M^{-1}(F_0 - C\dot{u}_0 - Ku_0)$
选择参数 $β\beta$ 和 $γ\gamma$

2. 计算等效刚度矩阵

$K^=K+γβΔtC+1βΔt2M\hat{K} = K + \frac{\gamma}{\beta\Delta t}C + \frac{1}{\beta\Delta t^2}M$

3. 时间步进循环

对于每个时间步 $n = 0, 1, 2, ...$ ：

a) 预测：

$u~n+1=un+Δt⋅u˙n+(12−β)Δt2⋅u¨n\tilde{u}_{n+1} = u_n + \Delta t \cdot \dot{u}_n + \left(\frac{1}{2}-\beta\right)\Delta t^2 \cdot \ddot{u}_n$

$u˙~n+1=u˙n+(1−γ)Δt⋅u¨n\tilde{\dot{u}}_{n+1} = \dot{u}_n + (1-\gamma)\Delta t \cdot \ddot{u}_n$

b) 计算等效荷载：

$F^n+1=Fn+1+M(1βΔt2u~n+1)+C(γβΔtu~n+1)\hat{F}_{n+1} = F_{n+1} + M\left(\frac{1}{\beta\Delta t^2}\tilde{u}_{n+1}\right) + C\left(\frac{\gamma}{\beta\Delta t}\tilde{u}_{n+1}\right)$

c) 求解位移：

$K^un+1=F^n+1\hat{K}u_{n+1} = \hat{F}_{n+1}$

d) 校正速度和加速度：

$u¨n+1=1βΔt2(un+1−u~n+1)\ddot{u}_{n+1} = \frac{1}{\beta\Delta t^2}(u_{n+1} - \tilde{u}_{n+1})$

$u˙n+1=u˙~n+1+γΔt⋅u¨n+1\dot{u}_{n+1} = \tilde{\dot{u}}_{n+1} + \gamma\Delta t \cdot \ddot{u}_{n+1}$

3.2.2 稳定性分析

Newmark-β法的稳定性取决于放大矩阵的谱半径。对于单自由度系统，放大矩阵为：

$\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}$

稳定性条件为谱半径 $ρ(A)≤1\rho(A) \leq 1$ 。

无条件稳定条件：

$γ≥0.5,β≥14(0.5+γ)2\gamma \geq 0.5, \quad \beta \geq \frac{1}{4}(0.5 + \gamma)^2$

条件稳定时的临界步长：

当 $γ=0.5\gamma = 0.5$ ， $β<0.25\beta < 0.25$ 时：

$Δtcr=1ωnγ/2−β\Delta t_{cr} = \frac{1}{\omega_n\sqrt{\gamma/2 - \beta}}$

3.3 其他隐式积分方法

3.3.1 Wilson-θ法

Wilson-θ法是一种无条件稳定的隐式方法，通过在时间区间 $t+\theta\Delta t]$ 内假设加速度线性变化来实现无条件稳定。

参数 $θ≥1.37\theta \geq 1.37$ （通常取1.4）保证无条件稳定。

3.3.2 Houbolt法

Houbolt法是一种多步法，利用前三步的位移来近似当前步的加速度和速度：

$u¨n+1=2un+1−5un+4un−1−un−2Δt2\ddot{u}_{n+1} = \frac{2u_{n+1} - 5u_n + 4u_{n-1} - u_{n-2}}{\Delta t^2}$

$u˙n+1=11un+1−18un+9un−1−2un−26Δt\dot{u}_{n+1} = \frac{11u_{n+1} - 18u_n + 9u_{n-1} - 2u_{n-2}}{6\Delta t}$

Houbolt法无条件稳定，但引入较大的数值阻尼。

4. 模态叠加法

4.1 模态叠加法基本原理

模态叠加法（Mode Superposition Method）是求解线性结构动力响应的高效方法，特别适用于长时间历程分析。

4.1.1 模态展开

对于 $n$ 自由度系统，位移向量可以表示为模态振型的线性组合：

$\sum_{i=1}^{n} \phi_i q_i(t) = \Phi q(t)$

其中， $ϕi\phi_i$ 是第 $i$ 阶模态振型， $q_i(t)$ 是模态坐标。

4.1.2 解耦方程

将模态展开代入运动方程：

$MΦq¨+CΦq˙+KΦq=F(t)M\Phi\ddot{q} + C\Phi\dot{q} + K\Phi q = F(t)$

左乘 $ΦT\Phi^T$ ：

$ΦTMΦq¨+ΦTCΦq˙+ΦTKΦq=ΦTF(t)\Phi^T M \Phi \ddot{q} + \Phi^T C \Phi \dot{q} + \Phi^T K \Phi q = \Phi^T F(t)$

利用模态的正交性：

$ϕiTMϕj={Mi,i=j0,i≠j\phi_i^T M \phi_j = \begin{cases} M_i, & i=j \\ 0, & i \neq j \end{cases}$

$ϕiTKϕj={Ki=ωi2Mi,i=j0,i≠j\phi_i^T K \phi_j = \begin{cases} K_i = \omega_i^2 M_i, & i=j \\ 0, & i \neq j \end{cases}$

得到解耦的模态方程：

$Miq¨i+Ciq˙i+Kiqi=Pi(t),i=1,2,...,nM_i \ddot{q}_i + C_i \dot{q}_i + K_i q_i = P_i(t), \quad i = 1, 2, ..., n$

其中， $Pi(t)=ϕiTF(t)P_i(t) = \phi_i^T F(t)$ 是模态力。

4.2 模态截断

4.2.1 截断准则

实际计算中，通常只需要考虑前 $r$ 阶模态（ $\ll n$ ）：

$\approx \sum_{i=1}^{r} \phi_i q_i(t)$

常用的截断准则包括：

频率截断：只考虑频率低于截断频率 $ωc\omega_c$ 的模态

累积有效质量：当累积有效质量达到总质量的90%以上时截断

模态贡献：根据模态参与因子的大小选择重要模态

4.2.2 模态参与因子

第 $i$ 阶模态的参与因子定义为：

$Γi=ϕiTM⋅1ϕiTMϕi=LiMi\Gamma_i = \frac{\phi_i^T M \cdot 1}{\phi_i^T M \phi_i} = \frac{L_i}{M_i}$

其中， $1$ 是全1向量， $L_i$ 是模态参与系数。

有效质量：

$Meff,i=(ϕiTM⋅1)2ϕiTMϕi=Γi2MiM_{eff,i} = \frac{(\phi_i^T M \cdot 1)^2}{\phi_i^T M \phi_i} = \Gamma_i^2 M_i$

4.3 模态叠加法算法流程

步骤1：特征值分析

求解广义特征值问题： $Kϕ=ω2MϕK\phi = \omega^2 M\phi$
获得固有频率 $ωi\omega_i$ 和模态振型 $ϕi\phi_i$

步骤2：模态选择

根据截断准则选择 $r$ 个主要模态
构造模态矩阵 $Φr=[ϕ1,ϕ2,...,ϕr]\Phi_r = [\phi_1, \phi_2, ..., \phi_r]$

步骤3：计算模态力

$Pr(t)=ΦrTF(t)P_r(t) = \Phi_r^T F(t)$

步骤4：求解模态响应

对每个模态求解： $q¨i+2ξiωiq˙i+ωi2qi=Pi(t)/Mi\ddot{q}_i + 2\xi_i\omega_i\dot{q}_i + \omega_i^2 q_i = P_i(t)/M_i$

步骤5：叠加模态响应

$\Phi_r q_r(t)$

5. 子结构方法

5.1 子结构方法概述

子结构方法（Substructure Method）将大型复杂结构分解为若干个子结构，分别分析后再组装，可以显著降低计算规模。

5.1.1 子结构划分

根据结构特点和计算需求，可以按以下方式划分子结构：

几何划分：按结构的几何特征划分
功能划分：按结构的功能部件划分
材料划分：按不同材料区域划分

5.1.2 子结构界面

子结构之间的连接界面分为：

位移界面：界面位移作为连接变量
力界面：界面力作为连接变量
混合界面：同时考虑位移和力

5.2 静力凝聚法

5.2.1 基本原理

静力凝聚法（Guyan Reduction）通过消去内部自由度，只保留界面自由度，实现模型降阶。

将自由度分为保留自由度 $u_m$ （界面）和凝聚自由度 $u_s$ （内部）：

$[KmmKmsKsmKss][umus]=[FmFs]\begin{bmatrix} K_{mm} & K_{ms} \\ K_{sm} & K_{ss} \end{bmatrix} \begin{bmatrix} u_m \\ u_s \end{bmatrix} = \begin{bmatrix} F_m \\ F_s \end{bmatrix}$

从第二个方程：

$K_{sm}u_m + K_{ss}u_s = F_s$

假设内部自由度无惯性效应（静力假设）：

$u_s = -K_{ss}^{-1}K_{sm}u_m + K_{ss}^{-1}F_s$

定义凝聚矩阵：

$\begin{bmatrix} I \\ -K_{ss}^{-1}K_{sm} \end{bmatrix}$

凝聚后的刚度矩阵和质量矩阵：

$Kˉ=TTKT=Kmm−KmsKss−1Ksm\bar{K} = T^T K T = K_{mm} - K_{ms}K_{ss}^{-1}K_{sm}$

$Mˉ=TTMT\bar{M} = T^T M T$

5.3 Craig-Bampton方法

5.3.1 方法概述

Craig-Bampton方法是一种动态子结构方法，综合考虑了界面自由度和固定界面模态，精度高于静力凝聚法。

5.3.2 坐标变换

Craig-Bampton方法的坐标变换为：

$[umus]=[I0ΨcΦn][umqn]\begin{bmatrix} u_m \\ u_s \end{bmatrix} = \begin{bmatrix} I & 0 \\ \Psi_c & \Phi_n \end{bmatrix} \begin{bmatrix} u_m \\ q_n \end{bmatrix}$

其中：

$u_m$ ：界面自由度
$u_s$ ：内部自由度
$Ψc=−Kss−1Ksm\Psi_c = -K_{ss}^{-1}K_{sm}$ ：约束模态（静力影响系数）
$Φn\Phi_n$ ：固定界面模态（内部自由振动模态）
$q_n$ ：模态坐标

5.3.3 减缩矩阵

变换矩阵：

$\begin{bmatrix} I & 0 \\ \Psi_c & \Phi_n \end{bmatrix}$

减缩后的刚度矩阵和质量矩阵：

$K~=TTKT=[K~mm00Λn]\tilde{K} = T^T K T = \begin{bmatrix} \tilde{K}_{mm} & 0 \\ 0 & \Lambda_n \end{bmatrix}$

$M~=TTMT=[M~mmM~mnM~nmI]\tilde{M} = T^T M T = \begin{bmatrix} \tilde{M}_{mm} & \tilde{M}_{mn} \\ \tilde{M}_{nm} & I \end{bmatrix}$

其中：

$K~mm=Kmm+KmsΨc\tilde{K}_{mm} = K_{mm} + K_{ms}\Psi_c$ （界面刚度）
$Λn=diag(ωi2)\Lambda_n = \text{diag}(\omega_i^2)$ （固定界面频率平方）
$M~mm=Mmm+MmsΨc+ΨcTMsm+ΨcTMssΨc\tilde{M}_{mm} = M_{mm} + M_{ms}\Psi_c + \Psi_c^T M_{sm} + \Psi_c^T M_{ss}\Psi_c$
$M~mn=(Mms+ΨcTMss)Φn\tilde{M}_{mn} = (M_{ms} + \Psi_c^T M_{ss})\Phi_n$

6. Python仿真实现

6.1 案例1：有限差分法求解振动方程

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

plt.switch_backend('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================================
# 1. 中心差分法
# ============================================================================

def central_difference(m, c, k, F, u0, v0, dt, t_end):
    """
    中心差分法求解单自由度系统
    
    参数:
        m: 质量
        c: 阻尼系数
        k: 刚度
        F: 外力函数
        u0: 初始位移
        v0: 初始速度
        dt: 时间步长
        t_end: 终止时间
    
    返回:
        t: 时间数组
        u: 位移数组
        v: 速度数组
        a: 加速度数组
    """
    # 时间离散
    t = np.arange(0, t_end + dt, dt)
    n_steps = len(t)
    
    # 初始化
    u = np.zeros(n_steps)
    v = np.zeros(n_steps)
    a = np.zeros(n_steps)
    
    # 初始条件
    u[0] = u0
    a[0] = (F(0) - c*v0 - k*u0) / m
    
    # 计算u_{-1}（虚拟步）
    u_minus1 = u0 - dt*v0 + (dt**2/2)*a[0]
    
    # 系数
    a0 = m/dt**2 + c/(2*dt)
    a1 = k - 2*m/dt**2
    a2 = m/dt**2 - c/(2*dt)
    
    # 第一步特殊处理
    u[1] = (F(dt) - a1*u0 - a2*u_minus1) / a0
    
    # 时间步进
    for n in range(1, n_steps-1):
        u[n+1] = (F(t[n+1]) - a1*u[n] - a2*u[n-1]) / a0
    
    # 计算速度和加速度（中心差分）
    v[1:-1] = (u[2:] - u[:-2]) / (2*dt)
    a[1:-1] = (u[2:] - 2*u[1:-1] + u[:-2]) / dt**2
    
    # 边界处理
    v[0] = v0
    v[-1] = (u[-1] - u[-2]) / dt
    a[0] = a[0]
    a[-1] = (u[-1] - 2*u[-2] + u[-3]) / dt**2
    
    return t, u, v, a

def case1a_central_difference():
    """
    案例1a：中心差分法求解简谐荷载响应
    """
    print("=" * 70)
    print("案例1a：中心差分法求解简谐荷载响应")
    print("=" * 70)
    
    # 系统参数
    m = 1.0      # kg
    k = 100.0    # N/m
    xi = 0.05    # 阻尼比
    omega_n = np.sqrt(k/m)  # 固有频率
    c = 2 * xi * omega_n * m  # 阻尼系数
    
    # 简谐荷载
    omega_f = 0.8 * omega_n  # 激励频率
    F0 = 10.0  # 幅值
    F = lambda t: F0 * np.sin(omega_f * t)
    
    # 初始条件
    u0 = 0.0
    v0 = 0.0
    
    # 时间参数
    dt = 0.01  # 时间步长
    t_end = 10.0  # 终止时间
    
    print(f"\n系统参数:")
    print(f"  质量 m = {m} kg")
    print(f"  刚度 k = {k} N/m")
    print(f"  阻尼比 ξ = {xi}")
    print(f"  固有频率 ωn = {omega_n:.3f} rad/s")
    print(f"  激励频率 ωf = {omega_f:.3f} rad/s")
    print(f"  时间步长 dt = {dt} s")
    
    # 稳定性检查
    dt_critical = 2 / omega_n
    print(f"\n稳定性分析:")
    print(f"  临界步长 dt_cr = 2/ωn = {dt_critical:.4f} s")
    print(f"  当前步长 dt = {dt} s")
    print(f"  稳定性: {'满足' if dt < dt_critical else '不满足'}")
    
    # 中心差分求解
    t_cd, u_cd, v_cd, a_cd = central_difference(m, c, k, F, u0, v0, dt, t_end)
    
    # 精确解（稳态响应）
    H = 1 / np.sqrt((1 - (omega_f/omega_n)**2)**2 + (2*xi*omega_f/omega_n)**2)
    phi = np.arctan2(2*xi*omega_f/omega_n, 1 - (omega_f/omega_n)**2)
    u_exact = (F0/k) * H * np.sin(omega_f * t_cd - phi)
    
    # 绘制结果
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 位移时程
    ax = axes[0, 0]
    ax.plot(t_cd, u_cd, 'b-', linewidth=1.5, label='中心差分解')
    ax.plot(t_cd, u_exact, 'r--', linewidth=1.5, label='精确解（稳态）')
    ax.set_xlabel('时间 (s)', fontsize=11)
    ax.set_ylabel('位移 (m)', fontsize=11)
    ax.set_title('位移响应对比', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 速度时程
    ax = axes[0, 1]
    ax.plot(t_cd, v_cd, 'g-', linewidth=1.5)
    ax.set_xlabel('时间 (s)', fontsize=11)
    ax.set_ylabel('速度 (m/s)', fontsize=11)
    ax.set_title('速度响应', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 加速度时程
    ax = axes[1, 0]
    ax.plot(t_cd, a_cd, 'm-', linewidth=1.5)
    ax.set_xlabel('时间 (s)', fontsize=11)
    ax.set_ylabel('加速度 (m/s²)', fontsize=11)
    ax.set_title('加速度响应', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 误差分析
    ax = axes[1, 1]
    error = np.abs(u_cd - u_exact)
    ax.semilogy(t_cd, error, 'r-', linewidth=1.5)
    ax.set_xlabel('时间 (s)', fontsize=11)
    ax.set_ylabel('绝对误差 (m)', fontsize=11)
    ax.set_title('数值误差', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case1a_central_difference.png', dpi=150, bbox_inches='tight')
    print("\n图形已保存: case1a_central_difference.png")
    plt.close()
    
    # 计算最大误差
    max_error = np.max(error)
    print(f"\n精度分析:")
    print(f"  最大绝对误差: {max_error:.6e} m")
    print(f"  相对误差: {max_error/np.max(np.abs(u_exact))*100:.2f}%")

def case1b_step_size_effect():
    """
    案例1b：时间步长对精度和稳定性的影响
    """
    print("\n" + "=" * 70)
    print("案例1b：时间步长影响分析")
    print("=" * 70)
    
    # 系统参数
    m = 1.0
    k = 100.0
    xi = 0.05
    omega_n = np.sqrt(k/m)
    c = 2 * xi * omega_n * m
    
    # 简谐荷载
    omega_f = 0.5 * omega_n
    F0 = 10.0
    F = lambda t: F0 * np.sin(omega_f * t)
    
    # 初始条件
    u0 = 0.0
    v0 = 0.0
    t_end = 5.0
    
    # 不同时间步长
    dt_values = [0.001, 0.005, 0.01, 0.02, 0.05]
    dt_critical = 2 / omega_n
    
    print(f"\n临界步长: dt_cr = {dt_critical:.4f} s")
    print(f"\n{'步长 dt':<12} {'稳定性':<12} {'最大误差':<15} {'计算时间':<12}")
    print("-" * 55)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    colors = ['#E74C3C', '#3498DB', '#2ECC71', '#F39C12', '#9B59B6']
    
    for dt, color in zip(dt_values, colors):
        stable = "稳定" if dt < dt_critical else "不稳定"
        
        try:
            t, u, v, a = central_difference(m, c, k, F, u0, v0, dt, t_end)
            
            # 精确解
            H = 1 / np.sqrt((1 - (omega_f/omega_n)**2)**2 + (2*xi*omega_f/omega_n)**2)
            phi = np.arctan2(2*xi*omega_f/omega_n, 1 - (omega_f/omega_n)**2)
            u_exact = (F0/k) * H * np.sin(omega_f * t - phi)
            
            max_error = np.max(np.abs(u - u_exact))
            
            print(f"{dt:<12.4f} {stable:<12} {max_error:<15.6e} {len(t):<12}")
            
            axes[0].plot(t, u, color=color, linewidth=1.5, label=f'dt={dt}')
            
        except:
            print(f"{dt:<12.4f} {stable:<12} {'发散':<15} {'-':<12}")
    
    # 精确解
    t_fine = np.linspace(0, t_end, 1000)
    H = 1 / np.sqrt((1 - (omega_f/omega_n)**2)**2 + (2*xi*omega_f/omega_n)**2)
    phi = np.arctan2(2*xi*omega_f/omega_n, 1 - (omega_f/omega_n)**2)
    u_exact_fine = (F0/k) * H * np.sin(omega_f * t_fine - phi)
    axes[0].plot(t_fine, u_exact_fine, 'k--', linewidth=2, label='精确解')
    
    axes[0].set_xlabel('时间 (s)', fontsize=11)
    axes[0].set_ylabel('位移 (m)', fontsize=11)
    axes[0].set_title('不同步长的位移响应', fontsize=12, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 稳定性区域示意
    ax = axes[1]
    omega_dt = np.linspace(0.01, 3, 100)
    spectral_radius = np.abs(2 - omega_dt**2 - omega_dt*np.sqrt(omega_dt**2 - 4)) / 2
    spectral_radius[omega_dt <= 2] = 1.0  # 稳定区域
    
    ax.plot(omega_dt, spectral_radius, 'b-', linewidth=2)
    ax.axvline(x=2, color='r', linestyle='--', linewidth=2, label='稳定性边界')
    ax.fill_between([0, 2], [0, 0], [2, 2], alpha=0.2, color='green', label='稳定区域')
    ax.set_xlabel('ωΔt', fontsize=11)
    ax.set_ylabel('谱半径 ρ', fontsize=11)
    ax.set_title('中心差分法稳定性分析', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 2)
    
    plt.tight_layout()
    plt.savefig('case1b_step_size_effect.png', dpi=150, bbox_inches='tight')
    print("\n图形已保存: case1b_step_size_effect.png")
    plt.close()

# ============================================================================
# 2. 主程序
# ============================================================================

if __name__ == '__main__':
    print("\n" + "=" * 70)
    print("案例1：有限差分法求解振动方程")
    print("=" * 70)
    
    case1a_central_difference()
    case1b_step_size_effect()
    
    print("\n" + "=" * 70)
    print("案例1计算完成！")
    print("=" * 70)

6.2 案例2：Newmark积分法与稳定性分析

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv

plt.switch_backend('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================================
# 1. Newmark-beta法
# ============================================================================

def newmark_beta(m, c, k, F, u0, v0, dt, t_end, beta=0.25, gamma=0.5):
    """
    Newmark-beta法求解单自由度系统
    
    参数:
        m, c, k: 质量、阻尼、刚度
        F: 外力函数
        u0, v0: 初始位移和速度
        dt: 时间步长
        t_end: 终止时间
        beta, gamma: Newmark参数
    
    返回:
        t, u, v, a: 时间、位移、速度、加速度
    """
    # 时间离散
    t = np.arange(0, t_end + dt, dt)
    n_steps = len(t)
    
    # 初始化
    u = np.zeros(n_steps)
    v = np.zeros(n_steps)
    a = np.zeros(n_steps)
    
    # 初始条件
    u[0] = u0
    v[0] = v0
    a[0] = (F(0) - c*v0 - k*u0) / m
    
    # 等效刚度
    k_eff = k + gamma/(beta*dt)*c + 1/(beta*dt**2)*m
    
    # 时间步进
    for n in range(n_steps - 1):
        # 预测
        u_pred = u[n] + dt*v[n] + (0.5-beta)*dt**2*a[n]
        v_pred = v[n] + (1-gamma)*dt*a[n]
        
        # 等效荷载
        f_eff = F(t[n+1]) + m/(beta*dt**2)*u_pred + gamma*c/(beta*dt)*u_pred
        
        # 求解位移
        u[n+1] = f_eff / k_eff
        
        # 校正速度和加速度
        a[n+1] = (u[n+1] - u_pred) / (beta*dt**2)
        v[n+1] = v_pred + gamma*dt*a[n+1]
    
    return t, u, v, a

def case2a_newmark_comparison():
    """
    案例2a：不同Newmark参数的比较
    """
    print("=" * 70)
    print("案例2a：Newmark-beta法参数比较")
    print("=" * 70)
    
    # 系统参数
    m = 1.0
    k = 100.0
    xi = 0.05
    omega_n = np.sqrt(k/m)
    c = 2 * xi * omega_n * m
    
    # 简谐荷载
    omega_f = 0.8 * omega_n
    F0 = 10.0
    F = lambda t: F0 * np.sin(omega_f * t)
    
    # 初始条件
    u0 = 0.0
    v0 = 0.0
    
    # 时间参数
    dt = 0.05  # 较大步长以显示差异
    t_end = 10.0
    
    print(f"\n系统参数:")
    print(f"  固有频率 ωn = {omega_n:.3f} rad/s")
    print(f"  激励频率 ωf = {omega_f:.3f} rad/s")
    print(f"  时间步长 dt = {dt} s")
    
    # 不同Newmark参数
    methods = [
        ('平均加速度法', 0.25, 0.5),
        ('线性加速度法', 1/6, 0.5),
        ('Fox-Goodwin', 1/12, 0.5),
    ]
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    colors = ['#E74C3C', '#3498DB', '#2ECC71']
    
    print(f"\n{'方法':<20} {'β':<10} {'γ':<10} {'稳定性':<15}")
    print("-" * 60)
    
    for (name, beta, gamma), color in zip(methods, colors):
        # 稳定性判断
        if gamma >= 0.5 and beta >= 0.25*(0.5 + gamma)**2:
            stability = "无条件稳定"
        elif gamma >= 0.5:
            omega_dt_max = 1/np.sqrt(gamma/2 - beta) if gamma/2 > beta else float('inf')
            stability = f"条件稳定(ωΔt<{omega_dt_max:.2f})"
        else:
            stability = "不稳定"
        
        print(f"{name:<20} {beta:<10.4f} {gamma:<10.4f} {stability:<15}")
        
        t, u, v, a = newmark_beta(m, c, k, F, u0, v0, dt, t_end, beta, gamma)
        axes[0].plot(t, u, color=color, linewidth=1.5, label=f'{name}')
    
    # 精确解
    t_fine = np.linspace(0, t_end, 1000)
    H = 1 / np.sqrt((1 - (omega_f/omega_n)**2)**2 + (2*xi*omega_f/omega_n)**2)
    phi = np.arctan2(2*xi*omega_f/omega_n, 1 - (omega_f/omega_n)**2)
    u_exact = (F0/k) * H * np.sin(omega_f * t_fine - phi)
    axes[0].plot(t_fine, u_exact, 'k--', linewidth=2, label='精确解')
    
    axes[0].set_xlabel('时间 (s)', fontsize=11)
    axes[0].set_ylabel('位移 (m)', fontsize=11)
    axes[0].set_title('不同Newmark参数的响应', fontsize=12, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 稳定性区域
    ax = axes[1]
    gamma_vals = np.linspace(0.3, 0.8, 100)
    beta_unconditional = 0.25 * (0.5 + gamma_vals)**2
    
    ax.fill_between(gamma_vals, beta_unconditional, 0.6, alpha=0.3, color='green', label='无条件稳定')
    ax.fill_between(gamma_vals, 0, beta_unconditional, alpha=0.3, color='yellow', label='条件稳定')
    
    # 标记常用方法
    ax.plot(0.5, 0.25, 'ro', markersize=10, label='平均加速度法')
    ax.plot(0.5, 1/6, 'bs', markersize=10, label='线性加速度法')
    ax.plot(0.5, 1/12, 'g^', markersize=10, label='Fox-Goodwin')
    
    ax.set_xlabel('γ', fontsize=11)
    ax.set_ylabel('β', fontsize=11)
    ax.set_title('Newmark-beta法稳定性区域', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0.3, 0.8)
    ax.set_ylim(0, 0.6)
    
    plt.tight_layout()
    plt.savefig('case2a_newmark_comparison.png', dpi=150, bbox_inches='tight')
    print("\n图形已保存: case2a_newmark_comparison.png")
    plt.close()

def case2b_amplitude_decay():
    """
    案例2b：数值阻尼分析
    """
    print("\n" + "=" * 70)
    print("案例2b：数值阻尼分析")
    print("=" * 70)
    
    # 无阻尼系统
    m = 1.0
    k = 100.0
    c = 0.0  # 无阻尼
    omega_n = np.sqrt(k/m)
    
    # 初始条件（自由振动）
    u0 = 1.0
    v0 = 0.0
    F = lambda t: 0.0
    
    # 时间参数
    dt = 0.1
    t_end = 20.0
    
    print(f"\n系统参数:")
    print(f"  无阻尼自由振动")
    print(f"  固有频率 ωn = {omega_n:.3f} rad/s")
    print(f"  初始位移 u0 = {u0} m")
    
    # 不同gamma值（引入数值阻尼）
    gamma_values = [0.5, 0.55, 0.6]
    beta = 0.25
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    colors = ['#E74C3C', '#3498DB', '#2ECC71']
    
    print(f"\n{'γ':<10} {'数值阻尼':<15}")
    print("-" * 30)
    
    for gamma, color in zip(gamma_values, colors):
        t, u, v, a = newmark_beta(m, c, k, F, u0, v0, dt, t_end, beta, gamma)
        
        # 计算振幅衰减
        peaks_idx = []
        for i in range(1, len(u)-1):
            if u[i] > u[i-1] and u[i] > u[i+1]:
                peaks_idx.append(i)
        
        if len(peaks_idx) >= 2:
            amplitude_ratio = u[peaks_idx[-1]] / u[peaks_idx[0]]
            log_dec = -np.log(amplitude_ratio) / (len(peaks_idx) - 1)
            xi_num = log_dec / (2*np.pi)
        else:
            xi_num = 0.0
        
        print(f"{gamma:<10.2f} {xi_num:<15.4f}")
        
        axes[0].plot(t, u, color=color, linewidth=1.5, label=f'γ={gamma}')
    
    # 精确解（无阻尼）
    t_exact = np.linspace(0, t_end, 1000)
    u_exact = u0 * np.cos(omega_n * t_exact)
    axes[0].plot(t_exact, u_exact, 'k--', linewidth=2, label='精确解')
    
    axes[0].set_xlabel('时间 (s)', fontsize=11)
    axes[0].set_ylabel('位移 (m)', fontsize=11)
    axes[0].set_title('数值阻尼效应', fontsize=12, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 周期延长效应
    ax = axes[1]
    dt_values = np.linspace(0.01, 0.3, 50)
    period_ratio = []
    
    for dt_test in dt_values:
        omega_dt = omega_n * dt_test
        # Newmark法的周期延长近似公式
        if beta == 0.25 and gamma == 0.5:
            T_ratio = 1 + (1/12)*(omega_dt)**2
        else:
            T_ratio = 1.0
        period_ratio.append(T_ratio)
    
    ax.plot(dt_values * omega_n, period_ratio, 'b-', linewidth=2)
    ax.set_xlabel('ωnΔt', fontsize=11)
    ax.set_ylabel('T_num/T_exact', fontsize=11)
    ax.set_title('周期延长效应', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.axhline(y=1.0, color='r', linestyle='--', label='精确周期')
    ax.legend()
    
    plt.tight_layout()
    plt.savefig('case2b_amplitude_decay.png', dpi=150, bbox_inches='tight')
    print("\n图形已保存: case2b_amplitude_decay.png")
    plt.close()

# ============================================================================
# 2. 主程序
# ============================================================================

if __name__ == '__main__':
    print("\n" + "=" * 70)
    print("案例2：Newmark积分法与稳定性分析")
    print("=" * 70)
    
    case2a_newmark_comparison()
    case2b_amplitude_decay()
    
    print("\n" + "=" * 70)
    print("案例2计算完成！")
    print("=" * 70)

6.3 案例3：模态叠加法与振型截断

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

plt.switch_backend('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================================
# 1. 模态叠加法
# ============================================================================

def modal_analysis(M, K, n_modes=None):
    """
    模态分析
    
    参数:
        M: 质量矩阵
        K: 刚度矩阵
        n_modes: 计算的模态数（None表示全部）
    
    返回:
        omega: 固有频率
        phi: 模态振型（质量归一化）
    """
    if n_modes is None:
        n_modes = M.shape[0]
    
    # 求解广义特征值问题
    eigenvalues, eigenvectors = eigh(K, M, subset_by_index=[0, n_modes-1])
    
    # 固有频率
    omega = np.sqrt(eigenvalues)
    
    # 质量归一化
    phi = eigenvectors.copy()
    for i in range(n_modes):
        mass_norm = np.sqrt(phi[:, i].T @ M @ phi[:, i])
        phi[:, i] = phi[:, i] / mass_norm
    
    return omega, phi

def modal_superposition(M, K, F_func, phi, omega, xi, u0, v0, t):
    """
    模态叠加法求解响应
    
    参数:
        M, K: 质量、刚度矩阵
        F_func: 外力函数（返回向量）
        phi: 模态振型矩阵
        omega: 固有频率数组
        xi: 阻尼比（标量或数组）
        u0, v0: 初始位移和速度
        t: 时间数组
    
    返回:
        u: 位移响应矩阵（n_dof x n_time）
    """
    n_dof = M.shape[0]
    n_modes = len(omega)
    n_time = len(t)
    
    # 初始条件的模态坐标
    q0 = phi.T @ M @ u0
    qdot0 = phi.T @ M @ v0
    
    # 响应矩阵
    u = np.zeros((n_dof, n_time))
    
    # 对每个模态求解
    for i in range(n_modes):
        omega_i = omega[i]
        xi_i = xi if np.isscalar(xi) else xi[i]
        
        # 模态质量（质量归一化后为1）
        M_i = 1.0
        
        # 模态力
        P_i = np.array([phi[:, i] @ F_func(ti) for ti in t])
        
        # 求解单自由度响应（杜哈梅尔积分或数值积分）
        q_i = solve_sdof(omega_i, xi_i, P_i/M_i, q0[i], qdot0[i], t)
        
        # 叠加
        u += np.outer(phi[:, i], q_i)
    
    return u

def solve_sdof(omega_n, xi, P, q0, qdot0, t):
    """
    求解单自由度系统响应（简谐激励假设）
    """
    dt = t[1] - t[0]
    n = len(t)
    q = np.zeros(n)
    
    # 简化为Newmark法求解
    omega_d = omega_n * np.sqrt(1 - xi**2)
    
    # 初始条件
    q[0] = q0
    if n > 1:
        q[1] = q0 + dt*qdot0 + (dt**2/2)*(-2*xi*omega_n*qdot0 - omega_n**2*q0 + P[0])
    
    # 中心差分
    for i in range(1, n-1):
        a_i = P[i] - 2*xi*omega_n*(q[i]-q[i-1])/dt - omega_n**2*q[i]
        q[i+1] = 2*q[i] - q[i-1] + dt**2 * a_i
    
    return q

def case3a_modal_analysis():
    """
    案例3a：多自由度系统模态分析
    """
    print("=" * 70)
    print("案例3a：多自由度系统模态分析")
    print("=" * 70)
    
    # 3自由度系统（剪切型框架）
    m = 1.0  # 质量
    k = 100.0  # 刚度
    
    # 质量矩阵
    M = m * np.array([
        [1.0, 0, 0],
        [0, 1.0, 0],
        [0, 0, 1.0]
    ])
    
    # 刚度矩阵（三层剪切框架）
    K = k * np.array([
        [2, -1, 0],
        [-1, 2, -1],
        [0, -1, 1]
    ])
    
    print(f"\n系统参数:")
    print(f"  质量 m = {m}")
    print(f"  刚度 k = {k}")
    
    # 模态分析
    omega, phi = modal_analysis(M, K)
    
    print(f"\n固有频率:")
    for i in range(len(omega)):
        print(f"  第{i+1}阶: ω{i+1} = {omega[i]:.4f} rad/s ({omega[i]/(2*np.pi):.4f} Hz)")
    
    print(f"\n模态振型（质量归一化）:")
    for i in range(len(omega)):
        print(f"  第{i+1}阶: {phi[:, i]}")
    
    # 绘制模态
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 模态振型图
    ax = axes[0]
    floors = np.array([1, 2, 3])
    colors = ['#E74C3C', '#3498DB', '#2ECC71']
    
    for i in range(len(omega)):
        mode_shape = np.concatenate([[0], phi[:, i]])  # 地面为0
        floors_plot = np.concatenate([[0], floors])
        ax.plot(mode_shape, floors_plot, 'o-', color=colors[i], linewidth=2, 
                label=f'第{i+1}阶 (ω={omega[i]:.2f})')
    
    ax.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
    ax.set_xlabel('归一化位移', fontsize=11)
    ax.set_ylabel('楼层', fontsize=11)
    ax.set_title('模态振型', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 3.5)
    
    # 频率分布
    ax = axes[1]
    mode_nums = np.arange(1, len(omega)+1)
    ax.bar(mode_nums, omega, color=colors, alpha=0.7)
    ax.set_xlabel('模态阶数', fontsize=11)
    ax.set_ylabel('固有频率 (rad/s)', fontsize=11)
    ax.set_title('固有频率分布', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('case3a_modal_analysis.png', dpi=150, bbox_inches='tight')
    print("\n图形已保存: case3a_modal_analysis.png")
    plt.close()
    
    return M, K, omega, phi

def case3b_modal_truncation():
    """
    案例3b：模态截断误差分析
    """
    print("\n" + "=" * 70)
    print("案例3b：模态截断误差分析")
    print("=" * 70)
    
    # 10自由度系统
    n_dof = 10
    m = 1.0
    k = 100.0
    
    # 质量矩阵
    M = m * np.eye(n_dof)
    
    # 刚度矩阵（三对角）
    K = np.zeros((n_dof, n_dof))
    for i in range(n_dof):
        K[i, i] = 2 * k
        if i > 0:
            K[i, i-1] = -k
            K[i-1, i] = -k
    K[-1, -1] = k  # 顶层
    
    # 模态分析（全部模态）
    omega_full, phi_full = modal_analysis(M, K)
    
    print(f"\n系统参数:")
    print(f"  自由度 n = {n_dof}")
    print(f"  前5阶频率: {omega_full[:5]}")
    
    # 简谐荷载（作用于顶层）
    omega_f = omega_full[0] * 0.5  # 低于第一阶频率
    F0 = 1.0
    F_func = lambda t: np.array([0]*(n_dof-1) + [F0*np.sin(omega_f*t)])
    
    # 初始条件
    u0 = np.zeros(n_dof)
    v0 = np.zeros(n_dof)
    
    # 时间
    t = np.linspace(0, 10, 1000)
    
    # 不同截断模态数
    n_modes_list = [1, 2, 3, 5, 10]
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    colors = ['#E74C3C', '#3498DB', '#2ECC71', '#F39C12', '#9B59B6']
    
    # 精确解（全部模态）
    u_exact = modal_superposition(M, K, F_func, phi_full, omega_full, 0.05, u0, v0, t)
    
    print(f"\n{'模态数':<10} {'顶层最大位移':<15} {'相对误差':<15}")
    print("-" * 45)
    
    for n_modes, color in zip(n_modes_list, colors):
        # 截断模态
        omega_trunc = omega_full[:n_modes]
        phi_trunc = phi_full[:, :n_modes]
        
        # 模态叠加
        u_trunc = modal_superposition(M, K, F_func, phi_trunc, omega_trunc, 0.05, u0, v0, t)
        
        # 计算误差
        max_disp_exact = np.max(np.abs(u_exact[-1, :]))
        max_disp_trunc = np.max(np.abs(u_trunc[-1, :]))
        error = abs(max_disp_trunc - max_disp_exact) / max_disp_exact * 100
        
        print(f"{n_modes:<10} {max_disp_trunc:<15.6f} {error:<15.2f}%")
        
        axes[0].plot(t, u_trunc[-1, :], color=color, linewidth=1.5, label=f'{n_modes}阶模态')
    
    axes[0].plot(t, u_exact[-1, :], 'k--', linewidth=2, label='精确解')
    axes[0].set_xlabel('时间 (s)', fontsize=11)
    axes[0].set_ylabel('顶层位移 (m)', fontsize=11)
    axes[0].set_title('不同截断模态数的响应', fontsize=12, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 收敛曲线
    ax = axes[1]
    errors = []
    for n_modes in range(1, n_dof+1):
        omega_trunc = omega_full[:n_modes]
        phi_trunc = phi_full[:, :n_modes]
        u_trunc = modal_superposition(M, K, F_func, phi_trunc, omega_trunc, 0.05, u0, v0, t)
        
        rmse = np.sqrt(np.mean((u_trunc[-1, :] - u_exact[-1, :])**2))
        errors.append(rmse)
    
    ax.semilogy(range(1, n_dof+1), errors, 'b-o', linewidth=2, markersize=6)
    ax.set_xlabel('截断模态数', fontsize=11)
    ax.set_ylabel('RMSE', fontsize=11)
    ax.set_title('模态截断收敛曲线', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.axhline(y=1e-6, color='r', linestyle='--', label='容差')
    ax.legend()
    
    plt.tight_layout()
    plt.savefig('case3b_modal_truncation.png', dpi=150, bbox_inches='tight')
    print("\n图形已保存: case3b_modal_truncation.png")
    plt.close()

# ============================================================================
# 2. 主程序
# ============================================================================

if __name__ == '__main__':
    print("\n" + "=" * 70)
    print("案例3：模态叠加法与振型截断")
    print("=" * 70)
    
    case3a_modal_analysis()
    case3b_modal_truncation()
    
    print("\n" + "=" * 70)
    print("案例3计算完成！")
    print("=" * 70)

6.4 案例4：子结构方法与模型降阶

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh, inv

plt.switch_backend('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================================
# 1. 子结构方法
# ============================================================================

def guyan_reduction(K, M, master_dofs):
    """
    Guyan静力凝聚法
    
    参数:
        K: 刚度矩阵
        M: 质量矩阵
        master_dofs: 保留自由度索引
    
    返回:
        K_reduced: 凝聚刚度矩阵
        M_reduced: 凝聚质量矩阵
        T: 变换矩阵
    """
    n = K.shape[0]
    slave_dofs = [i for i in range(n) if i not in master_dofs]
    
    # 分块矩阵
    K_mm = K[np.ix_(master_dofs, master_dofs)]
    K_ms = K[np.ix_(master_dofs, slave_dofs)]
    K_sm = K[np.ix_(slave_dofs, master_dofs)]
    K_ss = K[np.ix_(slave_dofs, slave_dofs)]
    
    M_mm = M[np.ix_(master_dofs, master_dofs)]
    M_ms = M[np.ix_(master_dofs, slave_dofs)]
    M_sm = M[np.ix_(slave_dofs, master_dofs)]
    M_ss = M[np.ix_(slave_dofs, slave_dofs)]
    
    # 凝聚矩阵
    K_ss_inv = inv(K_ss)
    Psi_c = -K_ss_inv @ K_sm
    
    # 变换矩阵
    T = np.zeros((n, len(master_dofs)))
    T[master_dofs, :] = np.eye(len(master_dofs))
    T[slave_dofs, :] = Psi_c
    
    # 凝聚后的矩阵
    K_reduced = T.T @ K @ T
    M_reduced = T.T @ M @ T
    
    return K_reduced, M_reduced, T

def craig_bampton(K, M, master_dofs, n_modes):
    """
    Craig-Bampton动态子结构法
    
    参数:
        K: 刚度矩阵
        M: 质量矩阵
        master_dofs: 界面自由度索引
        n_modes: 固定界面模态数
    
    返回:
        K_cb: CB刚度矩阵
        M_cb: CB质量矩阵
        T_cb: CB变换矩阵
    """
    n = K.shape[0]
    slave_dofs = [i for i in range(n) if i not in master_dofs]
    
    # 分块矩阵
    K_mm = K[np.ix_(master_dofs, master_dofs)]
    K_ms = K[np.ix_(master_dofs, slave_dofs)]
    K_sm = K[np.ix_(slave_dofs, master_dofs)]
    K_ss = K[np.ix_(slave_dofs, slave_dofs)]
    
---