主题068：结构动力学不确定性量化

摘要

不确定性量化（Uncertainty Quantification, UQ）是现代结构动力学分析的重要组成部分。本主题系统介绍结构动力学中不确定性的来源、分类和量化方法，包括随机有限元方法、多项式混沌展开、贝叶斯模型更新和可靠性分析等核心技术。通过Python实现多个典型案例，展示如何在实际工程问题中处理材料参数、几何尺寸、边界条件和载荷的不确定性，为结构的安全评估和优化设计提供理论基础和实践指导。

关键词

不确定性量化、随机有限元、多项式混沌、贝叶斯推断、可靠性分析、蒙特卡洛方法、结构动力学

---

1. 引言

1.1 不确定性的来源

在结构动力学问题中，不确定性主要来源于以下几个方面：

1. 材料参数的不确定性
材料属性如弹性模量、密度、阻尼比等存在固有的随机性。即使是同一批材料，其力学性能也会有一定的分散性。这种不确定性通常可以用概率分布来描述，如正态分布、对数正态分布等。

2. 几何尺寸的不确定性
制造公差、施工误差导致实际结构的几何尺寸与设计值存在偏差。对于复杂结构，这种几何不确定性可能显著影响动力响应特性。

3. 边界条件的不确定性
实际结构的支撑条件往往难以精确确定，如地基刚度、连接刚度等都可能存在不确定性。

4. 载荷的不确定性
环境载荷如风载、地震载荷、波浪载荷等具有显著的随机性和时空变异性。

5. 模型不确定性
数学模型对物理现实的简化、模型参数的选择、数值方法的近似等都会引入模型不确定性。

1.2 不确定性量化的意义

不确定性量化在结构动力学中的重要性体现在：

• 安全评估：考虑不确定性后，可以更准确地评估结构的失效概率
• 优化设计：在不确定性约束下进行鲁棒性优化设计
• 模型验证：通过试验数据更新模型参数，提高预测精度
• 决策支持：为工程决策提供量化的风险信息

1.3 不确定性分类

根据认知水平，不确定性可分为：

1. 偶然不确定性（Aleatory Uncertainty）
源于物理系统的固有随机性，不可通过增加数据消除。如材料性能的微观差异、环境载荷的随机波动等。

2. 认知不确定性（Epistemic Uncertainty）
源于知识的不完备，可通过增加数据或改进模型减少。如模型形式的选择、参数估计的误差等。

---

2. 概率论基础

2.1 随机变量与概率分布

定义：随机变量是定义在样本空间上的实值函数，用于描述随机现象。

概率密度函数（PDF）：
连续随机变量 $X$ 的概率密度函数 $f_X(x)$ 满足：

$$P(a\le X\le b)={\int }_a^b{f}_X(x)dx$$

累积分布函数（CDF）：
$$F_X(x)=P(X\le x)={\int }_{-\infty }^x{f}_X(t)dt$$

2.2 常用概率分布

正态分布（高斯分布）：
$$f_X(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

其中 $\mu$ 为均值，$\sigma$ 为标准差。

对数正态分布：
若 $\ln X$ 服从正态分布，则 $X$ 服从对数正态分布：
$$f_X(x)=\frac{1}{x\sigma \sqrt{2\pi }}\exp \left(-\frac{(\ln x-\mu {)}^2}{2{\sigma }^2}\right),x>0$$

对数正态分布常用于描述非负参数，如材料强度、疲劳寿命等。

均匀分布：
$$f_X(x)=\{\begin{array}{ll}\frac{1}{b-a},& a\le x\le b\\ 0,& \text{其他}\end{array}$$

当仅知道参数的上下界而缺乏更多信息时，常采用均匀分布。

威布尔分布：
$$f_X(x)=\frac{k}{\lambda }{\left(\frac{x}{\lambda }\right)}^{k-1}\exp \left(-{\left(\frac{x}{\lambda }\right)}^k\right),x\ge 0$$

威布尔分布广泛用于可靠性分析和寿命预测。

2.3 随机变量的数字特征

期望值（均值）：
$$E[X]={\int }_{-\infty }^{\infty }x{f}_X(x)dx$$

方差：
$$\text{Var}(X)=E[(X-E[X]{)}^2]={\int }_{-\infty }^{\infty }(x-E[X]{)}^2{f}_X(x)dx$$

标准差：
$${\sigma }_X=\sqrt{\text{Var}(X)}$$

变异系数：
$${\delta }_X=\frac{{\sigma }_X}{E[X]}$$

变异系数是无量纲的离散程度度量，便于比较不同量级参数的分散性。

---

3. 随机有限元方法

3.1 方法概述

随机有限元方法（Stochastic Finite Element Method, SFEM）是传统有限元方法的扩展，用于处理具有随机参数的结构分析问题。

基本思想：将随机参数和响应表示为随机变量的函数，通过统计方法求解响应的概率特性。

3.2 摄动法

摄动法是SFEM中最常用的方法之一，通过泰勒展开近似处理随机性。

一阶摄动法：
设结构参数 $a$ 为随机变量，可表示为：
$$a=\bar{a}+\Delta a$$

其中 $\bar{a}$ 为均值，$\Delta a$ 为随机扰动。

响应 $u$ 也可展开为：
$$u=\bar{u}+\Delta u$$

代入控制方程并保留一阶项，可得：
$$K(\bar{a})\bar{u}=F$$
$$K(\bar{a})\Delta u=-\frac{\partial K}{\partial a}\Delta a\cdot \bar{u}$$

响应统计特性：
$$E[u]\approx \bar{u}$$
$$\text{Var}(u)\approx E[\Delta u^2]$$

3.3 Neumann展开法

Neumann展开法是另一种处理随机刚度矩阵的方法。

将随机刚度矩阵表示为：
$$K=\bar{K}+\Delta K$$

其中 $\bar{K}$ 为均值刚度矩阵，$\Delta K$ 为随机扰动。

利用Neumann级数展开：
$$K^{-1}=(\bar{K}+\Delta K{)}^{-1}={\bar{K}}^{-1}-{\bar{K}}^{-1}\Delta K{\bar{K}}^{-1}+{\bar{K}}^{-1}\Delta K{\bar{K}}^{-1}\Delta K{\bar{K}}^{-1}-\cdots$$

响应的统计特性可通过级数截断近似计算。

3.4 谱随机有限元方法

谱随机有限元方法（Spectral SFEM）基于多项式混沌展开，将随机响应表示为随机变量的多项式函数。

Karhunen-Loève展开：
对于随机场 $H(x,\theta )$，可展开为：
$$H(x,\theta )=\bar{H}(x)+\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\varphi }_i(x){\xi }_i(\theta )$$

其中 ${\lambda }_i$ 和 ${\varphi }_i(x)$ 是协方差函数的特征值和特征函数，${\xi }_i$ 为标准正态随机变量。

---

4. 多项式混沌展开

4.1 混沌多项式基础

多项式混沌展开（Polynomial Chaos Expansion, PCE）由Wiener于1938年提出，Ghanem和Spanos将其应用于随机有限元分析。

基本思想：将随机响应表示为标准正态随机变量的多项式函数：
$$Y(\xi )=\sum _{i=0}^P{y}_i{\Psi }_i(\xi )$$

其中 ${\Psi }_i$ 为正交多项式（混沌多项式），$y_i$ 为确定性系数。

4.2 正交多项式族

不同类型的随机变量对应不同的正交多项式族：

随机变量类型	正交多项式	支撑集
高斯分布	Hermite多项式	$(-\infty ,\infty )$
伽马分布	Laguerre多项式	$[0,\infty )$
贝塔分布	Jacobi多项式	$[a,b]$
均匀分布	Legendre多项式	$[a,b]$

Hermite多项式：
$$H_n(\xi )=(-1{)}^n{e}^{{\xi }^2/2}\frac{d^n}{d{\xi }^n}{e}^{-{\xi }^2/2}$$

前几项：

• $H_0(\xi )=1$
• $H_1(\xi )=\xi$
• $H_2(\xi )={\xi }^2-1$
• $H_3(\xi )={\xi }^3-3\xi$

4.3 系数求解方法

1. 投影法（Galerkin投影）：
$$y_i=\frac{E[Y{\Psi }_i]}{E[{\Psi }_i^2]}$$

2. 配点法（Collocation）：
在选定的配点上求解确定性问题，通过插值得到PCE系数。

3. 回归法：
利用最小二乘法拟合PCE系数。

4.4 统计特性提取

一旦获得PCE系数，响应的统计特性可直接计算：

均值：
$$E[Y]=y_0$$

方差：
$$\text{Var}(Y)=\sum _{i=1}^P{y}_i^{2}E[{\Psi }_i^2]$$

高阶矩：
可通过多项式展开的高阶统计计算。

---

5. 贝叶斯模型更新

5.1 贝叶斯定理

贝叶斯定理是贝叶斯推断的基础：
$$P(\theta |D)=\frac{P(D|\theta )P(\theta )}{P(D)}$$

其中：

• $P(\theta )$：先验分布，表示观测数据前的参数信念
• $P(D|\theta )$：似然函数，表示在给定参数下观测数据的概率
• $P(\theta |D)$：后验分布，表示观测数据后的参数信念
• $P(D)$：证据，归一化常数

5.2 先验分布的选择

无信息先验：当缺乏先验知识时，采用均匀分布或Jeffreys先验。

信息先验：利用专家知识或历史数据构造先验分布。

共轭先验：选择与似然函数共轭的先验，使后验分布具有相同形式，便于解析计算。

5.3 似然函数的构造

假设观测误差服从正态分布：
$$D=Y(\theta )+ϵ,ϵ\sim N(0,{\sigma }^2)$$

似然函数为：
$$P(D|\theta )=\prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{(D_i-Y_i(\theta ){)}^2}{2{\sigma }^2}\right)$$

5.4 马尔可夫链蒙特卡洛（MCMC）

当后验分布复杂时，采用MCMC方法进行采样。

Metropolis-Hastings算法：

1. 初始化 ${\theta }^{(0)}$
2. 对于 $t=1,2,\cdots ,N$：
   • 从提议分布 $q({\theta }^{\ast }|{\theta }^{(t-1)})$ 采样候选点
   • 计算接受概率：$\alpha =\min \left(1,\frac{P({\theta }^{\ast }|D)q({\theta }^{(t-1)}|{\theta }^{\ast })}{P({\theta }^{(t-1)}|D)q({\theta }^{\ast }|{\theta }^{(t-1)})}\right)$
   • 以概率 $\alpha$ 接受 ${\theta }^{\ast }$，否则保持 ${\theta }^{(t-1)}$

收敛诊断：

• Gelman-Rubin统计量
• Geweke诊断
• 迹线图分析

---

6. 可靠性分析

6.1 失效概率的定义

结构的失效概率定义为：
$$P_f=P(g(X)\le 0)={\int }_{g(x)\le 0}{f}_X(x)dx$$

其中 $g(X)$ 为极限状态函数，$g(X)>0$ 表示安全，$g(X)\le 0$ 表示失效。

6.2 一次二阶矩法（FORM）

FORM通过线性近似计算失效概率。

标准正态变换：
将相关非正态随机变量 $X$ 变换为标准正态变量 $U$：
$$U=T(X)$$

可靠度指标：
在标准正态空间中，可靠度指标 $\beta$ 为原点到极限状态曲面的最短距离：
$$\beta =\underset{g(u)=0}{\min }\Vert u\Vert$$

失效概率近似：
$$P_f\approx \Phi (-\beta )$$

其中 $\Phi$ 为标准正态累积分布函数。

6.3 二次二阶矩法（SORM）

SORM在FORM的基础上，采用二次曲面近似极限状态函数，提高计算精度。

$$P_f\approx \Phi (-\beta )\prod _{i=1}^{n-1}(1-\beta {\kappa }_i{)}^{-1/2}$$

其中 ${\kappa }_i$ 为极限状态曲面在设计点处的主曲率。

6.4 重要性抽样

重要性抽样通过改变抽样分布，提高失效概率估计的效率。

$$P_f=\int I(g(x)\le 0)f_X(x)dx=\int I(g(x)\le 0)\frac{f_X(x)}{h(x)}h(x)dx$$

其中 $h(x)$ 为重要性抽样密度，通常选择以设计点为中心的正态分布。

---

7. Python实现案例

7.1 案例1：随机有限元分析

问题描述：
考虑一个悬臂梁，其弹性模量和截面惯性矩为随机变量。使用摄动法分析梁端位移的统计特性。

理论分析：
悬臂梁端部静挠度：
$$\delta =\frac{P{L}^3}{3EI}$$

当 $E$ 和 $I$ 为随机变量时，通过一阶摄动法：
$$\delta \approx \bar{\delta }+\frac{\partial \delta }{\partial E}\Delta E+\frac{\partial \delta }{\partial I}\Delta I$$

Python实现：
见 case1_stochastic_fem.py

7.2 案例2：多项式混沌展开

问题描述：
使用多项式混沌展开近似单自由度系统的随机响应，对比不同阶数PCE的精度。

理论分析：
SDOF系统运动方程：
$$m\ddot{u}+c\dot{u}+ku=F(t)$$

当 $k$ 和 $c$ 为随机变量时，响应 $u$ 可展开为：
$$u(t,\xi )=\sum _{i=0}^P{u}_i(t){\Psi }_i(\xi )$$

Python实现：
见 case2_polynomial_chaos.py

7.3 案例3：贝叶斯模型更新

问题描述：
基于观测的结构动力响应数据，使用MCMC方法更新结构参数的后验分布。

理论分析：
假设观测数据为结构频率，模型参数为刚度和质量。通过贝叶斯定理更新参数信念。

Python实现：
见 case3_bayesian_update.py

7.4 案例4：可靠性分析

问题描述：
使用FORM和蒙特卡洛方法计算结构在随机载荷作用下的失效概率。

理论分析：
定义极限状态函数：
$$g(R,S)=R-S$$

其中 $R$ 为结构抗力，$S$ 为载荷效应。

Python实现：
见 case4_reliability_analysis.py

---

8. 结果分析与讨论

8.1 随机有限元结果

通过摄动法和蒙特卡洛模拟的对比，验证摄动法在小幅不确定性情况下的有效性。当变异系数较大时，需要采用高阶摄动或蒙特卡洛方法。

8.2 多项式混沌展开精度

PCE的精度随多项式阶数增加而提高，但计算成本也随之增加。对于工程应用，通常2-4阶PCE即可满足精度要求。

8.3 贝叶斯更新的收敛性

MCMC方法的收敛速度受提议分布、参数相关性和后验分布复杂度的影响。通过调整提议分布的协方差矩阵，可以提高采样效率。

8.4 可靠性分析精度

FORM适用于线性或弱非线性极限状态函数，SORM适用于中等非线性情况。对于高度非线性问题，建议采用蒙特卡洛或重要性抽样方法。

---

9. 工程应用与展望

9.1 工程应用

不确定性量化方法已广泛应用于：

• 航空航天结构的可靠性设计
• 土木工程结构的安全评估
• 汽车结构的耐撞性分析
• 海上平台的疲劳寿命预测

9.2 发展趋势

1. 高维不确定性量化：
针对高维随机输入，发展稀疏网格、降维技术等方法。

2. 多尺度不确定性传递：
从材料微观到结构宏观的不确定性传递分析。

3. 数字孪生与实时更新：
结合监测数据，实现结构状态的实时不确定性更新。

4. 机器学习融合：
利用深度学习加速不确定性传播和可靠性分析。

---

10. 总结

本主题系统介绍了结构动力学不确定性量化的基本理论和方法，包括：

1. 随机有限元方法：摄动法、Neumann展开、谱方法
2. 多项式混沌展开：正交多项式、系数求解、统计特性提取
3. 贝叶斯模型更新：先验选择、MCMC采样、后验分析
4. 可靠性分析：FORM/SORM、重要性抽样、失效概率计算

通过Python实现四个典型案例，展示了不确定性量化方法在实际工程问题中的应用。这些方法为考虑不确定性的结构分析和设计提供了强有力的工具，有助于提高工程结构的安全性和经济性。

---

---

完整Python代码实现

以下是本主题的完整Python仿真代码：

import matplotlib
matplotlib.use('Agg')
"""
案例1：随机有限元分析
使用摄动法和蒙特卡洛方法分析悬臂梁在随机参数下的响应
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib.animation import FuncAnimation
import time
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 使用Agg后端
plt.switch_backend('Agg')


def cantilever_deflection(P, L, E, I):
    """
    计算悬臂梁端部挠度
    delta = P*L^3 / (3*E*I)
    """
    return P * L**3 / (3 * E * I)


def perturbation_analysis(P, L, E_mean, E_std, I_mean, I_std):
    """
    一阶摄动法分析
    """
    # 均值响应
    delta_mean = cantilever_deflection(P, L, E_mean, I_mean)
    
    # 偏导数
    d_delta_dE = -P * L**3 / (3 * E_mean**2 * I_mean)
    d_delta_dI = -P * L**3 / (3 * E_mean * I_mean**2)
    
    # 响应方差（假设E和I独立）
    var_delta = (d_delta_dE * E_std)**2 + (d_delta_dI * I_std)**2
    std_delta = np.sqrt(var_delta)
    
    return {
        'mean': delta_mean,
        'std': std_delta,
        'var': var_delta,
        'd_dE': d_delta_dE,
        'd_dI': d_delta_dI
    }


def monte_carlo_analysis(P, L, E_mean, E_std, I_mean, I_std, n_samples=10000):
    """
    蒙特卡洛模拟
    """
    # 生成随机样本
    E_samples = np.random.normal(E_mean, E_std, n_samples)
    I_samples = np.random.normal(I_mean, I_std, n_samples)
    
    # 确保物理合理性
    E_samples = np.clip(E_samples, E_mean * 0.5, E_mean * 2.0)
    I_samples = np.clip(I_samples, I_mean * 0.5, I_mean * 2.0)
    
    # 计算响应
    delta_samples = cantilever_deflection(P, L, E_samples, I_samples)
    
    return {
        'samples': delta_samples,
        'mean': np.mean(delta_samples),
        'std': np.std(delta_samples),
        'min': np.min(delta_samples),
        'max': np.max(delta_samples),
        'p5': np.percentile(delta_samples, 5),
        'p95': np.percentile(delta_samples, 95)
    }


def plot_comparison(perturbation_results, mc_results, P, L, E_mean, E_std, I_mean, I_std):
    """绘制摄动法和蒙特卡洛结果对比"""
    print("\n绘制摄动法与蒙特卡洛对比...")
    
    fig = plt.figure(figsize=(16, 10))
    gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)
    
    # 1. 统计特性对比
    ax1 = fig.add_subplot(gs[0, 0])
    methods = ['Perturbation', 'Monte Carlo']
    means = [perturbation_results['mean'], mc_results['mean']]
    stds = [perturbation_results['std'], mc_results['std']]
    
    x = np.arange(len(methods))
    width = 0.35
    
    bars1 = ax1.bar(x - width/2, means, width, label='Mean', color='steelblue', alpha=0.8)
    ax1_twin = ax1.twinx()
    bars2 = ax1_twin.bar(x + width/2, stds, width, label='Std', color='coral', alpha=0.8)
    
    ax1.set_ylabel('Mean Deflection (m)', color='steelblue', fontsize=10)
    ax1_twin.set_ylabel('Std Deflection (m)', color='coral', fontsize=10)
    ax1.set_xticks(x)
    ax1.set_xticklabels(methods)
    ax1.set_title('Statistical Comparison', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 添加数值标签
    for bar, val in zip(bars1, means):
        ax1.text(bar.get_x() + bar.get_width()/2., bar.get_height(),
                f'{val:.4f}', ha='center', va='bottom', fontsize=9)
    for bar, val in zip(bars2, stds):
        ax1_twin.text(bar.get_x() + bar.get_width()/2., bar.get_height(),
                     f'{val:.4f}', ha='center', va='bottom', fontsize=9)
    
    # 2. 蒙特卡洛样本分布
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.hist(mc_results['samples'], bins=50, density=True, alpha=0.7, 
            color='steelblue', edgecolor='black')
    
    # 拟合正态分布
    mu, std = stats.norm.fit(mc_results['samples'])
    x_fit = np.linspace(mc_results['samples'].min(), mc_results['samples'].max(), 100)
    ax2.plot(x_fit, stats.norm.pdf(x_fit, mu, std), 'r-', linewidth=2, label='Normal Fit')
    
    # 标记统计量
    ax2.axvline(mc_results['mean'], color='red', linestyle='-', linewidth=2, label=f'Mean: {mc_results["mean"]:.4f}m')
    ax2.axvline(mc_results['p5'], color='green', linestyle='--', linewidth=2, label=f'5%: {mc_results["p5"]:.4f}m')
    ax2.axvline(mc_results['p95'], color='orange', linestyle='--', linewidth=2, label=f'95%: {mc_results["p95"]:.4f}m')
    
    ax2.set_xlabel('Deflection (m)', fontsize=10)
    ax2.set_ylabel('Probability Density', fontsize=10)
    ax2.set_title('MC Distribution', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=8)
    ax2.grid(True, alpha=0.3)
    
    # 3. 参数敏感性分析
    ax3 = fig.add_subplot(gs[0, 2])
    params = ['E', 'I']
    sensitivities = [abs(perturbation_results['d_dE']) * E_std / perturbation_results['std'],
                    abs(perturbation_results['d_dI']) * I_std / perturbation_results['std']]
    
    colors = ['steelblue', 'coral']
    bars = ax3.bar(params, sensitivities, color=colors, alpha=0.8, edgecolor='black')
    ax3.set_ylabel('Sensitivity Index', fontsize=10)
    ax3.set_title('Parameter Sensitivity', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3, axis='y')
    
    for bar, val in zip(bars, sensitivities):
        ax3.text(bar.get_x() + bar.get_width()/2., bar.get_height(),
                f'{val:.2f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
    
    # 4. 收敛性分析
    ax4 = fig.add_subplot(gs[1, 0])
    sample_sizes = [100, 500, 1000, 2000, 5000, 10000]
    means_conv = []
    stds_conv = []
    
    for n in sample_sizes:
        subset = mc_results['samples'][:n]
        means_conv.append(np.mean(subset))
        stds_conv.append(np.std(subset))
    
    ax4.plot(sample_sizes, means_conv, 'b-o', linewidth=2, markersize=6, label='Mean')
    ax4.axhline(y=perturbation_results['mean'], color='r', linestyle='--', 
                linewidth=2, label=f'Perturbation Mean')
    ax4.set_xlabel('Sample Size', fontsize=10)
    ax4.set_ylabel('Mean Deflection (m)', fontsize=10)
    ax4.set_title('Convergence of Mean', fontsize=12, fontweight='bold')
    ax4.legend(fontsize=9)
    ax4.grid(True, alpha=0.3)
    ax4.set_xscale('log')
    
    # 5. 方差收敛
    ax5 = fig.add_subplot(gs[1, 1])
    ax5.plot(sample_sizes, stds_conv, 'g-s', linewidth=2, markersize=6, label='MC Std')
    ax5.axhline(y=perturbation_results['std'], color='r', linestyle='--', 
                linewidth=2, label=f'Perturbation Std')
    ax5.set_xlabel('Sample Size', fontsize=10)
    ax5.set_ylabel('Std Deflection (m)', fontsize=10)
    ax5.set_title('Convergence of Std', fontsize=12, fontweight='bold')
    ax5.legend(fontsize=9)
    ax5.grid(True, alpha=0.3)
    ax5.set_xscale('log')
    
    # 6. 概率密度对比
    ax6 = fig.add_subplot(gs[1, 2])
    
    # 摄动法假设的正态分布
    x_pert = np.linspace(perturbation_results['mean'] - 3*perturbation_results['std'],
                         perturbation_results['mean'] + 3*perturbation_results['std'], 100)
    y_pert = stats.norm.pdf(x_pert, perturbation_results['mean'], perturbation_results['std'])
    
    # 蒙特卡洛的核密度估计
    from scipy.stats import gaussian_kde
    kde = gaussian_kde(mc_results['samples'])
    x_mc = np.linspace(mc_results['samples'].min(), mc_results['samples'].max(), 100)
    y_mc = kde(x_mc)
    
    ax6.plot(x_pert, y_pert, 'b-', linewidth=2, label='Perturbation (Normal)')
    ax6.plot(x_mc, y_mc, 'r--', linewidth=2, label='MC (KDE)')
    ax6.set_xlabel('Deflection (m)', fontsize=10)
    ax6.set_ylabel('PDF', fontsize=10)
    ax6.set_title('PDF Comparison', fontsize=12, fontweight='bold')
    ax6.legend(fontsize=9)
    ax6.grid(True, alpha=0.3)
    
    plt.savefig('stochastic_fem_comparison.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("  Saved: stochastic_fem_comparison.png")


def plot_parameter_variation(P, L, E_mean, E_std, I_mean, I_std):
    """绘制参数变化对响应的影响"""
    print("\n绘制参数变化影响...")
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # E的变化影响
    ax1 = axes[0]
    E_range = np.linspace(E_mean - 3*E_std, E_mean + 3*E_std, 100)
    delta_E = cantilever_deflection(P, L, E_range, I_mean)
    
    ax1.plot(E_range/1e9, delta_E*1000, 'b-', linewidth=2)
    ax1.axvline(E_mean/1e9, color='r', linestyle='--', linewidth=2, label=f'Mean E={E_mean/1e9:.1f}GPa')
    ax1.axvline((E_mean-E_std)/1e9, color='orange', linestyle=':', linewidth=2, label=f'±1σ')
    ax1.axvline((E_mean+E_std)/1e9, color='orange', linestyle=':', linewidth=2)
    ax1.fill_between(E_range/1e9, delta_E*1000, alpha=0.3, color='lightblue')
    
    ax1.set_xlabel('Elastic Modulus E (GPa)', fontsize=11)
    ax1.set_ylabel('Deflection (mm)', fontsize=11)
    ax1.set_title('Effect of E Variation', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    
    # I的变化影响
    ax2 = axes[1]
    I_range = np.linspace(I_mean - 3*I_std, I_mean + 3*I_std, 100)
    delta_I = cantilever_deflection(P, L, E_mean, I_range)
    
    ax2.plot(I_range*1e6, delta_I*1000, 'r-', linewidth=2)
    ax2.axvline(I_mean*1e6, color='b', linestyle='--', linewidth=2, label=f'Mean I={I_mean*1e6:.1f}e-6 m⁴')
    ax2.axvline((I_mean-I_std)*1e6, color='orange', linestyle=':', linewidth=2, label=f'±1σ')
    ax2.axvline((I_mean+I_std)*1e6, color='orange', linestyle=':', linewidth=2)
    ax2.fill_between(I_range*1e6, delta_I*1000, alpha=0.3, color='lightcoral')
    
    ax2.set_xlabel('Moment of Inertia I (10⁻⁶ m⁴)', fontsize=11)
    ax2.set_ylabel('Deflection (mm)', fontsize=11)
    ax2.set_title('Effect of I Variation', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('parameter_variation.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("  Saved: parameter_variation.png")


def create_uncertainty_animation(P, L, E_mean, E_std, I_mean, I_std):
    """创建不确定性传播动画"""
    print("\n创建不确定性传播动画...")
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 左图：参数空间
    ax1 = axes[0]
    ax1.set_xlim(E_mean/1e9 - 4*E_std/1e9, E_mean/1e9 + 4*E_std/1e9)
    ax1.set_ylim(I_mean*1e6 - 4*I_std*1e6, I_mean*1e6 + 4*I_std*1e6)
    ax1.set_xlabel('E (GPa)', fontsize=11)
    ax1.set_ylabel('I (10⁻⁶ m⁴)', fontsize=11)
    ax1.set_title('Parameter Space Sampling', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 绘制置信椭圆
    from matplotlib.patches import Ellipse
    ellipse = Ellipse((E_mean/1e9, I_mean*1e6), 2*E_std/1e9, 2*I_std*1e6, 
                      fill=False, color='red', linestyle='--', linewidth=2, label='1σ')
    ax1.add_patch(ellipse)
    ellipse2 = Ellipse((E_mean/1e9, I_mean*1e6), 4*E_std/1e9, 4*I_std*1e6, 
                       fill=False, color='orange', linestyle=':', linewidth=2, label='2σ')
    ax1.add_patch(ellipse2)
    ax1.legend(fontsize=9)
    
    scatter1 = ax1.scatter([], [], c=[], cmap='viridis', s=20, alpha=0.6)
    
    # 右图：响应分布
    ax2 = axes[1]
    ax2.set_xlim(0.01, 0.05)
    ax2.set_ylim(0, 100)
    ax2.set_xlabel('Deflection (m)', fontsize=11)
    ax2.set_ylabel('Frequency', fontsize=11)
    ax2.set_title('Response Distribution', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 生成样本
    np.random.seed(42)
    n_total = 1000
    E_samples = np.random.normal(E_mean, E_std, n_total)
    I_samples = np.random.normal(I_mean, I_std, n_total)
    delta_samples = cantilever_deflection(P, L, E_samples, I_samples)
    
    def init():
        scatter1.set_offsets(np.empty((0, 2)))
        ax2.clear()
        ax2.set_xlim(0.01, 0.05)
        ax2.set_ylim(0, 100)
        ax2.set_xlabel('Deflection (m)', fontsize=11)
        ax2.set_ylabel('Frequency', fontsize=11)
        ax2.set_title('Response Distribution', fontsize=12, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        return scatter1,
    
    def update(frame):
        n = (frame + 1) * 20
        
        # 更新参数空间
        x = E_samples[:n] / 1e9
        y = I_samples[:n] * 1e6
        c = delta_samples[:n] * 1000
        scatter1.set_offsets(np.c_[x, y])
        scatter1.set_array(c)
        
        # 更新响应分布
        ax2.clear()
        ax2.set_xlim(0.01, 0.05)
        ax2.set_ylim(0, 100)
        ax2.set_xlabel('Deflection (m)', fontsize=11)
        ax2.set_ylabel('Frequency', fontsize=11)
        ax2.set_title(f'Response Distribution (n={n})', fontsize=12, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        
        if n > 10:
            ax2.hist(delta_samples[:n], bins=30, alpha=0.7, color='steelblue', edgecolor='black')
        
        return scatter1,
    
    anim = FuncAnimation(fig, update, frames=50, init_func=init, blit=False, interval=100)
    anim.save('uncertainty_propagation.gif', writer='pillow', fps=10, dpi=100)
    plt.close()
    print("  Saved: uncertainty_propagation.gif")


def main():
    """主函数"""
    print("="*70)
    print("案例1：随机有限元分析")
    print("="*70)
    
    # 问题参数
    print("\n问题参数:")
    P = 10000.0  # 载荷 (N)
    L = 2.0      # 梁长度 (m)
    
    # 随机参数（弹性模量和惯性矩）
    E_mean = 2.0e11  # 弹性模量均值 (Pa)
    E_std = 0.1e11   # 弹性模量标准差 (Pa), 5%变异系数
    
    I_mean = 8.33e-6  # 惯性矩均值 (m⁴)
    I_std = 0.417e-6  # 惯性矩标准差 (m⁴), 5%变异系数
    
    print(f"  载荷 P: {P:.0f} N")
    print(f"  长度 L: {L:.1f} m")
    print(f"  弹性模量 E: {E_mean/1e9:.1f} ± {E_std/1e9:.1f} GPa")
    print(f"  惯性矩 I: {I_mean*1e6:.2f} ± {I_std*1e6:.2f} × 10⁻⁶ m⁴")
    
    # 确定性分析（均值参数）
    delta_det = cantilever_deflection(P, L, E_mean, I_mean)
    print(f"\n确定性分析结果:")
    print(f"  端部挠度: {delta_det*1000:.2f} mm")
    
    # 摄动法分析
    print("\n" + "="*50)
    print("摄动法分析")
    print("="*50)
    
    t_start = time.time()
    perturbation_results = perturbation_analysis(P, L, E_mean, E_std, I_mean, I_std)
    t_pert = time.time() - t_start
    
    print(f"  计算时间: {t_pert*1000:.2f} ms")
    print(f"  均值挠度: {perturbation_results['mean']*1000:.2f} mm")
    print(f"  标准差: {perturbation_results['std']*1000:.2f} mm")
    print(f"  变异系数: {perturbation_results['std']/perturbation_results['mean']*100:.1f}%")
    print(f"  对E的敏感度: {perturbation_results['d_dE']:.6f}")
    print(f"  对I的敏感度: {perturbation_results['d_dI']:.6f}")
    
    # 蒙特卡洛模拟
    print("\n" + "="*50)
    print("蒙特卡洛模拟")
    print("="*50)
    
    n_samples = 10000
    print(f"  样本数: {n_samples}")
    
    t_start = time.time()
    mc_results = monte_carlo_analysis(P, L, E_mean, E_std, I_mean, I_std, n_samples)
    t_mc = time.time() - t_start
    
    print(f"  计算时间: {t_mc:.2f} s")
    print(f"  均值挠度: {mc_results['mean']*1000:.2f} mm")
    print(f"  标准差: {mc_results['std']*1000:.2f} mm")
    print(f"  最小值: {mc_results['min']*1000:.2f} mm")
    print(f"  最大值: {mc_results['max']*1000:.2f} mm")
    print(f"  5%分位数: {mc_results['p5']*1000:.2f} mm")
    print(f"  95%分位数: {mc_results['p95']*1000:.2f} mm")
    
    # 方法对比
    print("\n" + "="*50)
    print("方法对比")
    print("="*50)
    
    mean_error = abs(perturbation_results['mean'] - mc_results['mean']) / mc_results['mean'] * 100
    std_error = abs(perturbation_results['std'] - mc_results['std']) / mc_results['std'] * 100
    speedup = t_mc / t_pert
    
    print(f"  均值误差: {mean_error:.2f}%")
    print(f"  标准差误差: {std_error:.2f}%")
    print(f"  计算加速比: {speedup:.0f}x")
    
    # 生成可视化
    print("\n" + "="*50)
    print("生成可视化结果")
    print("="*50)
    
    plot_comparison(perturbation_results, mc_results, P, L, E_mean, E_std, I_mean, I_std)
    plot_parameter_variation(P, L, E_mean, E_std, I_mean, I_std)
    create_uncertainty_animation(P, L, E_mean, E_std, I_mean, I_std)
    
    # 总结
    print("\n" + "="*70)
    print("案例1完成！")
    print("="*70)
    print("\n主要结论:")
    print("  1. 摄动法计算快速，适合小幅不确定性")
    print("  2. 蒙特卡洛方法精度高，但计算成本高")
    print("  3. 两种方法的结果吻合良好")
    print("  4. 参数变异系数5%导致响应变异系数约7%")
    print("\n生成的文件:")
    print("  - stochastic_fem_comparison.png: 方法对比")
    print("  - parameter_variation.png: 参数影响")
    print("  - uncertainty_propagation.gif: 不确定性传播动画")


if __name__ == "__main__":
    main()

代码说明

1. 参数设置区：定义系统的质量、刚度、阻尼等基本参数
2. 核心计算模块：实现振动方程的求解算法
3. 可视化模块：生成分析图表和动画
4. 结果输出：保存图片文件供文档使用

运行上述代码将生成本主题所需的所有可视化素材。