主题097: 仿真结果验证与确认(V&V)

目录

1. 概述
2. 验证与确认的基本概念
3. 代码验证方法
4. 模型确认方法
5. 不确定性量化在V&V中的应用
6. 工程案例分析

---

概述

为什么需要V&V？

在工程仿真领域，验证(Verification)与确认(Validation)是确保仿真结果可信度的核心过程。随着计算机仿真技术在航空航天、汽车、能源等关键领域的广泛应用，仿真结果的可靠性直接关系到产品安全、性能和成本。

V&V的重要性体现在：

1. 确保仿真可信度：通过系统化的验证和确认流程，量化仿真结果的不确定性
2. 降低工程风险：在物理试验前识别模型缺陷，避免昂贵的设计返工
3. 满足行业标准：符合ASME V&V 10、AIAA G-077等国际标准要求
4. 支持决策制定：为工程设计和管理决策提供可靠的数据支撑

V&V的发展历程

第一阶段(1960s-1980s)：主要关注数值方法的正确性验证

• 有限元方法的数学基础验证
• 数值算法稳定性分析
• 网格收敛性研究

第二阶段(1990s-2000s)：模型确认成为研究热点

• 物理试验与仿真对比方法
• 误差估计和不确定性量化
• ASME V&V 10标准发布(2006)

第三阶段(2010s-至今)：系统化V&V框架

• 数字孪生中的V&V应用
• 机器学习辅助的模型确认
• 实时验证与自适应模型更新

---

验证与确认的基本概念

定义与区别

验证(Verification)：

“我们是否正确地求解了数学模型？”

验证关注的是数值解与数学模型之间的一致性，回答"代码是否正确实现了数学模型"这一问题。

确认(Validation)：

“我们是否求解了正确的数学模型？”

确认关注的是数学模型与物理现实之间的一致性，回答"模型是否能够准确预测真实物理现象"这一问题。

V&V的数学框架

验证的数学表达：

$${ϵ}_{num}=||u_{exact}-u_h||$$

其中：

• $u_{exact}$：数学模型的精确解
• $u_h$：数值解（离散化解）
• ${ϵ}_{num}$：数值误差

确认的数学表达：

$${ϵ}_{model}=||y_{exp}-y_{sim}||$$

其中：

• $y_{exp}$：试验观测值
• $y_{sim}$：仿真预测值
• ${ϵ}_{model}$：模型误差

V&V流程图

┌─────────────────────────────────────────────────────────────┐
│                     真实物理系统                            │
└──────────────────────┬──────────────────────────────────────┘
                       │ 概念建模
                       ▼
┌─────────────────────────────────────────────────────────────┐
│                     数学模型                                │
│         (PDEs, ODEs, Integral Equations)                   │
└──────────────────────┬──────────────────────────────────────┘
                       │ 离散化
                       ▼
┌─────────────────────────────────────────────────────────────┐
│                   计算机模型                                │
│            (FEM, FDM, BEM等数值方法)                       │
└──────────────────────┬──────────────────────────────────────┘
                       │ 代码实现
                       ▼
┌─────────────────────────────────────────────────────────────┐
│                   仿真代码                                  │
└─────────────────────────────────────────────────────────────┘

验证路径：数学模型 → 离散化 → 代码实现 → 数值解
确认路径：真实物理系统 → 数学模型 → 仿真结果 → 试验对比

---

代码验证方法

1. 制造解方法(Method of Manufactured Solutions, MMS)

MMS是验证数值代码最严格的方法之一。其核心思想是：

1. 假设一个解析解 $u_{manufactured}$
2. 将其代入控制方程，计算源项
3. 用数值代码求解带源项的方程
4. 对比数值解与制造解

实施步骤：

# MMS示例：验证一维热传导求解器

# 步骤1: 选择制造解
u_exact = sin(pi*x) * exp(-t)

# 步骤2: 计算源项
# du/dt = d²u/dx² + f
# f = du/dt - d²u/dx²
f = -pi**2 * sin(pi*x) * exp(-t) - pi**2 * sin(pi*x) * exp(-t)

# 步骤3: 数值求解
u_numerical = solve_heat_equation(f, boundary_conditions)

# 步骤4: 计算误差
error = norm(u_numerical - u_exact)

MMS的优势：

• 可以验证代码的所有项（边界条件、源项、非线性项）
• 可以精确计算误差，进行收敛性分析
• 适用于复杂方程组

2. 网格收敛性分析

网格收敛性分析是验证数值方法精度的基本手段。

Richardson外推法：

假设数值解具有如下渐近展开：

$$u_h=u_{exact}+C{h}^p+O(h^{p+1})$$

其中：

• $h$：网格尺寸
• $p$：收敛阶数
• $C$：常数

网格收敛指标(GCI)：

$$GCI=\frac{F_s|ϵ|}{r^p-1}$$

其中：

• $F_s$：安全因子（通常取1.25）
• $ϵ=(u_{fine}-u_{coarse})/u_{fine}$：相对变化
• $r$：网格细化比
• $p$：理论收敛阶

实施流程：

1. 在三个不同网格密度下计算（粗、中、细）
2. 计算观测收敛阶数：

$$p_{obs}=\frac{\ln \left(\frac{u_{medium}-u_{coarse}}{u_{fine}-u_{medium}}\right)}{\ln (r)}$$

3. 验证 $p_{obs}\approx p_{theory}$
4. 计算GCI评估数值不确定性

3. 基准问题验证

使用公认的基准测试问题验证代码正确性。

结构动力学基准问题：

问题类型	基准案例	验证重点
模态分析	悬臂梁固有频率	特征值求解器
瞬态响应	简支梁冲击响应	时间积分算法
非线性	大变形梁弯曲	几何非线性处理
接触	赫兹接触问题	接触算法

---

模型确认方法

1. 确认实验设计

确认实验的基本原则：

1. 层次化方法：从简单到复杂逐步确认
2. 充分覆盖：实验条件覆盖仿真应用范围
3. 测量精度：试验不确定度小于仿真不确定度
4. 可重复性：实验结果具有统计意义

确认层次：

层次0: 单元级确认
    └── 材料本构关系验证

层次1: 子系统确认
    └── 组件级动力学特性

层次2: 系统级确认
    └── 整体结构响应

层次3: 全尺度确认
    └── 实际工况验证

2. 定量确认指标

均方根误差(RMSE)：

$$RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^N(y_{sim,i}-y_{exp,i}{)}^2}$$

决定系数(R²)：

$$R^2=1-\frac{\sum (y_{exp}-y_{sim}{)}^2}{\sum (y_{exp}-{\bar{y}}_{exp}{)}^2}$$

平均绝对百分比误差(MAPE)：

$$MAPE=\frac{100\%}{N}\sum _{i=1}^N\left|\frac{y_{exp,i}-y_{sim,i}}{y_{exp,i}}\right|$$

确认度量(Validation Metric)：

$$d=\frac{|{\bar{y}}_{sim}-{\bar{y}}_{exp}|}{\sqrt{u_{sim}^2+u_{exp}^2}}$$

其中 $u_{sim}$ 和 $u_{exp}$ 分别为仿真和试验的不确定度。

3. 模型修正与更新

贝叶斯模型修正：

$$P(\theta |D)=\frac{P(D|\theta )P(\theta )}{P(D)}$$

其中：

• $\theta$：模型参数
• $D$：试验数据
• $P(\theta |D)$：后验分布
• $P(D|\theta )$：似然函数
• $P(\theta )$：先验分布

卡尔曼滤波更新：

用于实时模型参数更新：

$${\theta }_{k+1}={\theta }_k+K_k(y_{exp,k}-y_{sim,k})$$

$$K_k=P_k{H}^T(H{P}_k{H}^T+R{)}^{-1}$$

---

不确定性量化在V&V中的应用

1. 误差来源分类

数值误差：

• 离散化误差（网格尺寸、时间步长）
• 迭代收敛误差
• 舍入误差

模型形式误差：

• 简化假设引入的误差
• 未建模物理现象
• 边界条件近似

参数不确定性：

• 材料参数分散性
• 几何尺寸公差
• 载荷不确定性

2. 总不确定度估计

总仿真不确定度：

$$u_{sim}^2=u_{num}^2+u_{model}^2+u_{param}^2$$

验证不确定度：

$$u_{ver}^2=u_{num}^2+u_{input}^2$$

确认不确定度：

$$u_{val}^2=u_{sim}^2+u_{exp}^2$$

3. 置信区间与预测区间

置信区间(CI)：反映参数估计的不确定性

$$CI=\hat{\theta }\pm t_{\alpha /2,n-1}\cdot SE(\hat{\theta })$$

预测区间(PI)：反映未来观测值的范围

$$PI=\hat{y}\pm t_{\alpha /2,n-1}\cdot \sqrt{MSE+SE(\hat{y}{)}^2}$$

---

工程案例分析

案例1: 代码验证与计算结果确认

目标：验证结构动力学求解器的正确性

方法：

1. MMS方法验证代码实现
2. 网格收敛性分析
3. 与解析解对比

预期成果：

• 验证数值方法的收敛阶数
• 量化数值误差
• 建立网格密度选择准则

import matplotlib
matplotlib.use('Agg')
"""
主题097 - 案例1: 代码验证与计算结果确认
使用制造解方法(MMS)和网格收敛性分析验证结构动力学求解器
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.patches import FancyBboxPatch, Rectangle, FancyArrowPatch

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


def analytical_solution_string_vibration(x, t, L=1.0, c=1.0, n_modes=5):
    """
    弦振动问题的解析解
    使用分离变量法得到的级数解
    """
    u = np.zeros_like(x)
    for n in range(1, n_modes + 1):
        # 初始条件: u(x,0) = sin(n*pi*x/L)
        # 解: u(x,t) = sin(n*pi*x/L) * cos(n*pi*c*t/L)
        An = 1.0 / n  # 振幅系数
        u += An * np.sin(n * np.pi * x / L) * np.cos(n * np.pi * c * t / L)
    return u


def finite_difference_solver(x, t, L=1.0, c=1.0):
    """
    有限差分法求解弦振动方程
    使用中心差分格式
    """
    dx = x[1] - x[0]
    dt = t[1] - t[0]
    r = c * dt / dx  # CFL数
    
    nx = len(x)
    nt = len(t)
    
    # 稳定性检查
    if r > 1:
        print(f"警告: CFL数 {r:.2f} > 1, 数值不稳定")
    
    # 初始化
    u = np.zeros((nt, nx))
    
    # 初始条件
    for i in range(nx):
        u[0, i] = np.sin(np.pi * x[i] / L)
    
    # 初始速度为0，使用特殊处理
    # u[1] = u[0] (因为初始速度为0)
    u[1, :] = u[0, :]
    
    # 时间推进
    for n in range(1, nt - 1):
        for i in range(1, nx - 1):
            u[n+1, i] = 2 * (1 - r**2) * u[n, i] + r**2 * (u[n, i+1] + u[n, i-1]) - u[n-1, i]
        
        # 边界条件 (固定端)
        u[n+1, 0] = 0
        u[n+1, -1] = 0
    
    return u


def manufactured_solution_verification():
    """
    制造解方法(MMS)验证
    验证数值代码是否正确实现了数学模型
    """
    print("=" * 70)
    print("案例1: 代码验证与计算结果确认")
    print("=" * 70)
    
    print("\n1. 制造解方法(MMS)验证")
    print("-" * 60)
    
    # 定义制造解: u = sin(pi*x) * cos(pi*t)
    def u_manufactured(x, t):
        return np.sin(np.pi * x) * np.cos(np.pi * t)
    
    # 计算源项（将制造解代入控制方程 u_tt = u_xx + f）
    # u_tt = -pi^2 * sin(pi*x) * cos(pi*t)
    # u_xx = -pi^2 * sin(pi*x) * cos(pi*t)
    # f = u_tt - u_xx = 0
    
    # 参数设置
    L = 1.0
    c = 1.0
    nx = 51
    nt = 101
    
    x = np.linspace(0, L, nx)
    t = np.linspace(0, 2.0, nt)
    
    print(f"网格设置: nx={nx}, nt={nt}")
    print(f"CFL数: {c * (t[1]-t[0]) / (x[1]-x[0]):.3f}")
    
    # 数值求解
    print("运行数值求解...")
    u_numerical = finite_difference_solver(x, t, L, c)
    
    # 计算误差
    print("\n计算误差...")
    errors = []
    max_errors = []
    
    for i, ti in enumerate(t):
        u_exact = u_manufactured(x, ti)
        error = np.abs(u_numerical[i] - u_exact)
        max_error = np.max(error)
        
        errors.append(error)
        max_errors.append(max_error)
    
    errors = np.array(errors)
    max_errors = np.array(max_errors)
    
    print(f"\nMMS验证结果:")
    print(f"  最大绝对误差: {np.max(max_errors):.6e}")
    print(f"  平均最大误差: {np.mean(max_errors):.6e}")
    print(f"  最终时刻误差: {max_errors[-1]:.6e}")
    
    return {
        'x': x,
        't': t,
        'u_numerical': u_numerical,
        'errors': max_errors,
        'error_dist': errors,
        'u_manufactured': u_manufactured
    }


def grid_convergence_analysis():
    """
    网格收敛性分析
    验证数值方法的收敛阶数
    """
    print("\n\n2. 网格收敛性分析")
    print("-" * 60)
    
    # 不同网格密度
    nx_list = [11, 21, 41, 81, 161]
    
    # 存储结果
    results = {
        'nx': [],
        'h': [],
        'max_errors': [],
        'l2_errors': []
    }
    
    L = 1.0
    c = 1.0
    t_final = 1.0
    
    print("\n网格收敛性分析结果:")
    print(f"{'网格数':<10} {'h':<12} {'最大误差':<15} {'L2误差':<15} {'误差比':<15}")
    print("-" * 70)
    
    prev_error = None
    
    for nx in nx_list:
        x = np.linspace(0, L, nx)
        dx = x[1] - x[0]
        
        # 根据CFL条件设置时间步长
        dt = 0.9 * dx / c
        nt = int(t_final / dt) + 1
        t = np.linspace(0, t_final, nt)
        
        # 数值求解
        u_numerical = finite_difference_solver(x, t, L, c)
        
        # 计算误差
        u_exact = np.sin(np.pi * x) * np.cos(np.pi * t[-1])
        error = np.abs(u_numerical[-1] - u_exact)
        max_error = np.max(error)
        l2_error = np.sqrt(np.mean(error**2))
        
        results['nx'].append(nx)
        results['h'].append(dx)
        results['max_errors'].append(max_error)
        results['l2_errors'].append(l2_error)
        
        # 计算误差比
        if prev_error is not None:
            ratio = prev_error / max_error
            print(f"{nx:<10} {dx:<12.4f} {max_error:<15.6e} {l2_error:<15.6e} {ratio:<15.2f}")
        else:
            print(f"{nx:<10} {dx:<12.4f} {max_error:<15.6e} {l2_error:<15.6e} {'-':<15}")
        
        prev_error = max_error
    
    # 计算观测收敛阶数
    print("\n观测收敛阶数:")
    h = np.array(results['h'])
    max_errors = np.array(results['max_errors'])
    
    for i in range(1, len(max_errors)):
        if max_errors[i] > 0 and max_errors[i-1] > 0:
            p_obs = np.log(max_errors[i-1] / max_errors[i]) / np.log(h[i-1] / h[i])
            print(f"  细化步骤 {i}: p_obs = {p_obs:.2f}")
    
    # 计算GCI
    print("\n网格收敛指标(GCI):")
    Fs = 1.25  # 安全因子
    r = 2.0    # 网格细化比
    p = 2.0    # 理论收敛阶
    
    for i in range(1, len(max_errors)):
        epsilon = abs(max_errors[i] - max_errors[i-1])
        GCI = Fs * epsilon / (r**p - 1)
        print(f"  步骤 {i}: GCI = {GCI:.6e}")
    
    return results


def visualize_mms_results(mms_data):
    """可视化MMS验证结果"""
    fig = plt.figure(figsize=(16, 5))
    
    x = mms_data['x']
    t = mms_data['t']
    u_num = mms_data['u_numerical']
    errors = mms_data['errors']
    error_dist = mms_data['error_dist']
    u_manu = mms_data['u_manufactured']
    
    # 子图1: 特定时刻的解对比
    ax1 = fig.add_subplot(1, 3, 1)
    time_indices = [0, 25, 50, 75, 100]
    colors = plt.cm.viridis(np.linspace(0, 1, len(time_indices)))
    
    for idx, color in zip(time_indices, colors):
        if idx < len(t):
            u_exact = u_manu(x, t[idx])
            ax1.plot(x, u_exact, '--', color=color, linewidth=2, 
                    label=f't={t[idx]:.2f}s (解析)', alpha=0.7)
            ax1.plot(x, u_num[idx], 'o', color=color, markersize=4,
                    label=f't={t[idx]:.2f}s (数值)', alpha=0.7)
    
    ax1.set_xlabel('位置 x (m)', fontsize=11)
    ax1.set_ylabel('位移 u (m)', fontsize=11)
    ax1.set_title('制造解 vs 数值解', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=8, ncol=2)
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 误差时程
    ax2 = fig.add_subplot(1, 3, 2)
    ax2.semilogy(t, errors, 'b-', linewidth=2)
    ax2.set_xlabel('时间 (s)', fontsize=11)
    ax2.set_ylabel('最大绝对误差', fontsize=11)
    ax2.set_title('数值误差演化', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    ax2.axhline(y=1e-3, color='r', linestyle='--', label='1e-3阈值')
    ax2.legend(fontsize=9)
    
    # 子图3: 空间误差分布
    ax3 = fig.add_subplot(1, 3, 3)
    final_error = error_dist[-1]
    ax3.bar(x, final_error, width=0.015, alpha=0.7, color='steelblue', edgecolor='black')
    ax3.set_xlabel('位置 x (m)', fontsize=11)
    ax3.set_ylabel('绝对误差 (m)', fontsize=11)
    ax3.set_title('最终时刻空间误差分布', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('案例1_MMS验证结果.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("\nMMS验证结果图已保存到: 案例1_MMS验证结果.png")
    plt.close()


def visualize_convergence(grid_data):
    """可视化网格收敛性分析结果"""
    fig = plt.figure(figsize=(16, 5))
    
    nx = np.array(grid_data['nx'])
    h = np.array(grid_data['h'])
    max_errors = np.array(grid_data['max_errors'])
    l2_errors = np.array(grid_data['l2_errors'])
    
    # 子图1: 误差随网格变化
    ax1 = fig.add_subplot(1, 3, 1)
    
    ax1.loglog(h, max_errors, 'o-', linewidth=2, markersize=8, label='最大误差')
    ax1.loglog(h, l2_errors, 's-', linewidth=2, markersize=8, label='L2误差')
    
    # 理论收敛线
    h_theory = np.linspace(h[-1], h[0], 100)
    ax1.loglog(h_theory, h_theory**2 * max_errors[-1] / h[-1]**2, 
              'k--', linewidth=2, label='理论2阶收敛', alpha=0.7)
    
    ax1.set_xlabel('网格尺寸 h (m)', fontsize=11)
    ax1.set_ylabel('误差', fontsize=11)
    ax1.set_title('网格收敛性分析', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3, which='both')
    
    # 子图2: 误差比
    ax2 = fig.add_subplot(1, 3, 2)
    
    error_ratios = []
    for i in range(1, len(max_errors)):
        ratio = max_errors[i-1] / max_errors[i]
        error_ratios.append(ratio)
    
    ax2.bar(range(1, len(error_ratios)+1), error_ratios, alpha=0.7, color='steelblue', edgecolor='black')
    ax2.axhline(y=4, color='r', linestyle='--', linewidth=2, label='理论值 (2^2=4)')
    ax2.set_xlabel('细化步骤', fontsize=11)
    ax2.set_ylabel('误差比', fontsize=11)
    ax2.set_title('误差收敛比', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 收敛阶数估计
    ax3 = fig.add_subplot(1, 3, 3)
    
    p_obs_list = []
    for i in range(1, len(max_errors)):
        if max_errors[i] > 0 and max_errors[i-1] > 0:
            p = np.log(max_errors[i-1] / max_errors[i]) / np.log(2)
            p_obs_list.append(p)
    
    ax3.plot(range(1, len(p_obs_list)+1), p_obs_list, 'o-', linewidth=2, markersize=10, color='steelblue')
    ax3.axhline(y=2, color='r', linestyle='--', linewidth=2, label='理论阶数 p=2', alpha=0.7)
    ax3.set_xlabel('细化步骤', fontsize=11)
    ax3.set_ylabel('观测收敛阶数', fontsize=11)
    ax3.set_title('收敛阶数验证', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    ax3.set_ylim(0, 4)
    
    plt.tight_layout()
    plt.savefig('案例1_网格收敛性分析.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("网格收敛性分析图已保存到: 案例1_网格收敛性分析.png")
    plt.close()


def create_vv_process_animation():
    """创建V&V流程动画"""
    print("\n正在生成V&V流程动画...")
    
    fig, ax = plt.subplots(figsize=(12, 8))
    ax.set_xlim(0, 10)
    ax.set_ylim(0, 10)
    ax.axis('off')
    ax.set_title('代码验证与确认流程', fontsize=16, fontweight='bold', pad=20)
    
    # 定义框的位置和内容
    boxes = [
        {'pos': (1, 8), 'size': (2.5, 1), 'text': '数学模型\n(PDEs)', 'color': 'lightblue'},
        {'pos': (4.5, 8), 'size': (2.5, 1), 'text': '离散化\n(FEM/FDM)', 'color': 'lightgreen'},
        {'pos': (7, 8), 'size': (2.5, 1), 'text': '数值代码\n(求解器)', 'color': 'lightyellow'},
        {'pos': (4.5, 5.5), 'size': (2.5, 1), 'text': '制造解\n(MMS)', 'color': 'lightcoral'},
        {'pos': (1, 3), 'size': (2.5, 1), 'text': '网格收敛\n分析', 'color': 'plum'},
        {'pos': (4.5, 3), 'size': (2.5, 1), 'text': '基准问题\n验证', 'color': 'lightsalmon'},
        {'pos': (7, 3), 'size': (2.5, 1), 'text': '误差估计\n与量化', 'color': 'lightgray'},
        {'pos': (4.5, 0.5), 'size': (2.5, 1), 'text': '验证完成\n代码可信', 'color': 'lightgreen'},
    ]
    
    # 绘制框
    box_patches = []
    text_objects = []
    
    for box in boxes:
        rect = FancyBboxPatch(box['pos'], box['size'][0], box['size'][1],
                             boxstyle="round,pad=0.1", 
                             facecolor=box['color'], 
                             edgecolor='black', 
                             linewidth=2,
                             alpha=0.3)
        ax.add_patch(rect)
        box_patches.append(rect)
        
        text = ax.text(box['pos'][0] + box['size'][0]/2, 
                      box['pos'][1] + box['size'][1]/2,
                      box['text'], 
                      ha='center', va='center', 
                      fontsize=11, fontweight='bold',
                      alpha=0.3)
        text_objects.append(text)
    
    # 添加箭头
    arrows = [
        ((2.25, 8), (4.5, 8)),      # 数学模型 -> 离散化
        ((5.25, 8), (7, 8)),        # 离散化 -> 数值代码
        ((7, 7.5), (5.75, 5.5)),    # 数值代码 -> MMS
        ((4.5, 5), (2.25, 3)),      # MMS -> 网格收敛
        ((5.75, 5), (5.75, 3)),     # MMS -> 基准问题
        ((7, 5.5), (8.25, 3)),      # MMS -> 误差估计
        ((2.25, 3), (4.5, 1.5)),    # 网格收敛 -> 验证完成
        ((5.75, 3), (5.75, 1.5)),   # 基准问题 -> 验证完成
        ((8.25, 3), (7, 1.5)),      # 误差估计 -> 验证完成
    ]
    
    arrow_objects = []
    for start, end in arrows:
        arrow = FancyArrowPatch(start, end, 
                               arrowstyle='->', 
                               mutation_scale=20,
                               linewidth=2, 
                               color='black',
                               alpha=0.3)
        ax.add_patch(arrow)
        arrow_objects.append(arrow)
    
    def init():
        return box_patches + text_objects + arrow_objects
    
    def update(frame):
        n_boxes = len(boxes)
        n_arrows = len(arrows)
        
        # 依次高亮框
        if frame < n_boxes:
            box_patches[frame].set_alpha(1)
            text_objects[frame].set_alpha(1)
        
        # 然后高亮箭头
        elif frame < n_boxes + n_arrows:
            arrow_idx = frame - n_boxes
            arrow_objects[arrow_idx].set_alpha(1)
        
        return box_patches + text_objects + arrow_objects
    
    total_frames = len(boxes) + len(arrows) + 5
    anim = FuncAnimation(fig, update, init_func=init, frames=total_frames,
                        interval=500, blit=False)
    
    writer = PillowWriter(fps=2)
    anim.save('案例1_VV流程动画.gif', writer=writer)
    print("动画已保存到: 案例1_VV流程动画.gif")
    plt.close()


def main():
    """主函数"""
    print("=" * 70)
    print("主题097 - 案例1: 代码验证与计算结果确认")
    print("=" * 70)
    
    # 1. MMS验证
    mms_data = manufactured_solution_verification()
    
    # 2. 网格收敛性分析
    grid_data = grid_convergence_analysis()
    
    # 生成可视化
    print("\n" + "="*60)
    print("生成可视化结果")
    print("="*60)
    
    visualize_mms_results(mms_data)
    visualize_convergence(grid_data)
    create_vv_process_animation()
    
    print("\n" + "=" * 70)
    print("案例1完成！所有结果已保存。")
    print("=" * 70)


if __name__ == "__main__":
    main()

案例2: 模型验证与不确定性量化

目标：评估仿真模型的预测能力

方法：

1. 参数敏感性分析
2. 不确定性传播分析
3. 模型可信度评估

预期成果：

• 识别关键影响参数
• 量化预测不确定度
• 提出模型改进建议

import matplotlib
matplotlib.use('Agg')
"""
主题097 - 案例2: 模型验证与不确定性量化
评估仿真模型的预测能力和参数敏感性
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.patches import FancyBboxPatch, Rectangle, Circle
from scipy import stats
from scipy.stats import norm, uniform

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


class BridgeModel:
    """简化的桥梁结构模型"""
    
    def __init__(self, L=100, E=30e9, I=10.0, rho=2500, A=20.0):
        """
        初始化桥梁模型
        
        参数:
            L: 跨度 (m)
            E: 弹性模量 (Pa)
            I: 截面惯性矩 (m^4)
            rho: 密度 (kg/m^3)
            A: 截面积 (m^2)
        """
        self.L = L
        self.E = E
        self.I = I
        self.rho = rho
        self.A = A
        self.m = rho * A  # 单位长度质量
    
    def natural_frequency(self, n):
        """计算第n阶固有频率 (Hz)"""
        # 简支梁固有频率公式
        return (n * np.pi)**2 / (2 * np.pi * self.L**2) * np.sqrt(self.E * self.I / self.m)
    
    def deflection_static(self, P, a):
        """计算集中载荷作用下的静挠度"""
        # 简支梁中点集中载荷挠度公式
        x = self.L / 2
        if x <= a:
            delta = P * x * (self.L - a) * (2 * a * self.L - a**2 - x**2) / (6 * self.E * self.I * self.L)
        else:
            delta = P * a * (self.L - x) * (2 * x * self.L - x**2 - a**2) / (6 * self.E * self.I * self.L)
        return delta
    
    def max_deflection(self, P):
        """中点集中载荷最大挠度"""
        return P * self.L**3 / (48 * self.E * self.I)
    
    def fundamental_period(self):
        """计算基本周期"""
        return 1 / self.natural_frequency(1)


def monte_carlo_simulation(model, param_distributions, n_samples=1000):
    """
    蒙特卡洛模拟进行不确定性传播分析
    
    参数:
        model: 结构模型
        param_distributions: 参数概率分布字典
        n_samples: 样本数量
    
    返回:
        仿真结果样本
    """
    print(f"\n运行蒙特卡洛模拟 (n={n_samples})...")
    
    # 存储结果
    frequencies = []
    deflections = []
    periods = []
    
    # 参数样本
    param_samples = {
        'E': [],
        'I': [],
        'rho': []
    }
    
    for i in range(n_samples):
        # 从分布中采样
        E_sample = param_distributions['E'].rvs()
        I_sample = param_distributions['I'].rvs()
        rho_sample = param_distributions['rho'].rvs()
        
        # 更新模型参数
        model.E = E_sample
        model.I = I_sample
        model.rho = rho_sample
        model.m = rho_sample * model.A
        
        # 计算响应
        freq = model.natural_frequency(1)
        deflection = model.max_deflection(10000)  # 10kN载荷
        period = model.fundamental_period()
        
        frequencies.append(freq)
        deflections.append(deflection)
        periods.append(period)
        
        param_samples['E'].append(E_sample)
        param_samples['I'].append(I_sample)
        param_samples['rho'].append(rho_sample)
        
        if (i + 1) % 200 == 0:
            print(f"  完成 {(i+1)/n_samples*100:.0f}%")
    
    return {
        'frequencies': np.array(frequencies),
        'deflections': np.array(deflections),
        'periods': np.array(periods),
        'param_samples': param_samples
    }


def sensitivity_analysis(model, param_ranges):
    """
    参数敏感性分析
    
    参数:
        model: 结构模型
        param_ranges: 参数变化范围字典
    
    返回:
        敏感性指标
    """
    print("\n进行参数敏感性分析...")
    
    # 基准值
    base_values = {
        'E': model.E,
        'I': model.I,
        'rho': model.rho
    }
    
    # 基准响应
    base_freq = model.natural_frequency(1)
    base_deflection = model.max_deflection(10000)
    
    sensitivities = {}
    
    for param_name, (min_val, max_val) in param_ranges.items():
        # 变化参数，保持其他参数不变
        values = np.linspace(min_val, max_val, 50)
        freq_responses = []
        deflection_responses = []
        
        for val in values:
            if param_name == 'E':
                model.E = val
            elif param_name == 'I':
                model.I = val
            elif param_name == 'rho':
                model.rho = val
                model.m = val * model.A
            
            freq_responses.append(model.natural_frequency(1))
            deflection_responses.append(model.max_deflection(10000))
        
        # 恢复基准值
        model.E = base_values['E']
        model.I = base_values['I']
        model.rho = base_values['rho']
        model.m = base_values['rho'] * model.A
        
        # 计算敏感性系数 (归一化)
        freq_sensitivity = np.std(freq_responses) / base_freq
        deflection_sensitivity = np.std(deflection_responses) / base_deflection
        
        sensitivities[param_name] = {
            'values': values,
            'frequencies': np.array(freq_responses),
            'deflections': np.array(deflection_responses),
            'freq_sensitivity': freq_sensitivity,
            'deflection_sensitivity': deflection_sensitivity
        }
    
    return sensitivities


def calculate_validation_metrics(sim_data, exp_data):
    """
    计算模型确认指标
    
    参数:
        sim_data: 仿真数据
        exp_data: 试验数据
    
    返回:
        确认指标字典
    """
    # 均方根误差
    rmse = np.sqrt(np.mean((sim_data - exp_data)**2))
    
    # 决定系数 R²
    ss_res = np.sum((exp_data - sim_data)**2)
    ss_tot = np.sum((exp_data - np.mean(exp_data))**2)
    r2 = 1 - ss_res / ss_tot
    
    # 平均绝对百分比误差
    mape = np.mean(np.abs((exp_data - sim_data) / exp_data)) * 100
    
    # 最大误差
    max_error = np.max(np.abs(sim_data - exp_data))
    
    return {
        'RMSE': rmse,
        'R2': r2,
        'MAPE': mape,
        'MaxError': max_error
    }


def main():
    """主函数"""
    print("=" * 70)
    print("主题097 - 案例2: 模型验证与不确定性量化")
    print("=" * 70)
    
    # 创建桥梁模型
    bridge = BridgeModel()
    
    print("\n桥梁模型参数:")
    print(f"  跨度 L = {bridge.L} m")
    print(f"  弹性模量 E = {bridge.E/1e9:.1f} GPa")
    print(f"  惯性矩 I = {bridge.I:.2f} m^4")
    print(f"  密度 rho = {bridge.rho} kg/m^3")
    
    # 基准响应
    print("\n基准响应:")
    print(f"  第1阶频率 = {bridge.natural_frequency(1):.4f} Hz")
    print(f"  第2阶频率 = {bridge.natural_frequency(2):.4f} Hz")
    print(f"  中点挠度(10kN) = {bridge.max_deflection(10000)*1000:.4f} mm")
    
    # 1. 不确定性量化
    print("\n" + "="*60)
    print("1. 不确定性量化")
    print("="*60)
    
    # 定义参数概率分布
    param_distributions = {
        'E': norm(loc=30e9, scale=3e9),      # 弹性模量: 均值30GPa, 标准差3GPa
        'I': norm(loc=10.0, scale=0.5),       # 惯性矩: 均值10, 标准差0.5
        'rho': norm(loc=2500, scale=100)      # 密度: 均值2500, 标准差100
    }
    
    # 蒙特卡洛模拟
    mc_results = monte_carlo_simulation(bridge, param_distributions, n_samples=1000)
    
    # 统计结果
    print("\n蒙特卡洛模拟统计结果:")
    print(f"{'响应量':<20} {'均值':<15} {'标准差':<15} {'变异系数':<15}")
    print("-" * 70)
    
    for name, data in [('基频(Hz)', mc_results['frequencies']), 
                       ('挠度(mm)', mc_results['deflections']*1000),
                       ('周期(s)', mc_results['periods'])]:
        mean_val = np.mean(data)
        std_val = np.std(data)
        cov = std_val / mean_val * 100
        print(f"{name:<20} {mean_val:<15.4f} {std_val:<15.4f} {cov:<15.2f}%")
    
    # 2. 敏感性分析
    print("\n" + "="*60)
    print("2. 参数敏感性分析")
    print("="*60)
    
    param_ranges = {
        'E': (24e9, 36e9),      # 弹性模量变化范围
        'I': (8.0, 12.0),        # 惯性矩变化范围
        'rho': (2200, 2800)      # 密度变化范围
    }
    
    sensitivities = sensitivity_analysis(bridge, param_ranges)
    
    print("\n参数敏感性指标 (变异系数):")
    print(f"{'参数':<10} {'频率敏感性':<15} {'挠度敏感性':<15}")
    print("-" * 45)
    for param, data in sensitivities.items():
        print(f"{param:<10} {data['freq_sensitivity']*100:<15.2f}% {data['deflection_sensitivity']*100:<15.2f}%")
    
    # 3. 模型确认
    print("\n" + "="*60)
    print("3. 模型确认指标")
    print("="*60)
    
    # 模拟试验数据（添加噪声）
    np.random.seed(42)
    n_exp = 20
    
    # 理论值
    freq_theory = bridge.natural_frequency(1)
    deflection_theory = bridge.max_deflection(10000)
    
    # 试验数据（添加5%随机误差）
    freq_exp = freq_theory * (1 + 0.05 * np.random.randn(n_exp))
    deflection_exp = deflection_theory * (1 + 0.05 * np.random.randn(n_exp))
    
    # 仿真数据（蒙特卡洛均值）
    freq_sim = np.mean(mc_results['frequencies'])
    deflection_sim = np.mean(mc_results['deflections'])
    
    # 计算确认指标
    freq_metrics = calculate_validation_metrics(
        np.full_like(freq_exp, freq_sim), 
        freq_exp
    )
    
    deflection_metrics = calculate_validation_metrics(
        np.full_like(deflection_exp, deflection_sim), 
        deflection_exp
    )
    
    print("\n频率确认指标:")
    for metric, value in freq_metrics.items():
        print(f"  {metric}: {value:.4f}")
    
    print("\n挠度确认指标:")
    for metric, value in deflection_metrics.items():
        print(f"  {metric}: {value:.4f}")
    
    # 生成可视化
    print("\n" + "="*60)
    print("生成可视化结果")
    print("="*60)
    
    visualize_uncertainty(mc_results)
    visualize_sensitivity(sensitivities)
    visualize_validation(mc_results, freq_exp, deflection_exp)
    create_uncertainty_animation(mc_results)
    
    print("\n" + "=" * 70)
    print("案例2完成！所有结果已保存。")
    print("=" * 70)


def visualize_uncertainty(mc_results):
    """可视化不确定性分析结果"""
    fig = plt.figure(figsize=(16, 5))
    
    # 子图1: 基频分布
    ax1 = fig.add_subplot(1, 3, 1)
    freqs = mc_results['frequencies']
    ax1.hist(freqs, bins=30, density=True, alpha=0.7, color='steelblue', edgecolor='black')
    
    # 拟合正态分布
    mu, std = norm.fit(freqs)
    x = np.linspace(freqs.min(), freqs.max(), 100)
    ax1.plot(x, norm.pdf(x, mu, std), 'r-', linewidth=2, label=f'正态拟合\nμ={mu:.3f}, σ={std:.3f}')
    
    ax1.axvline(mu, color='g', linestyle='--', linewidth=2, label=f'均值={mu:.3f}')
    ax1.axvline(mu - 2*std, color='orange', linestyle=':', linewidth=1.5, alpha=0.7)
    ax1.axvline(mu + 2*std, color='orange', linestyle=':', linewidth=1.5, alpha=0.7, label='95%置信区间')
    
    ax1.set_xlabel('基频 (Hz)', fontsize=11)
    ax1.set_ylabel('概率密度', fontsize=11)
    ax1.set_title('基频概率分布', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 挠度分布
    ax2 = fig.add_subplot(1, 3, 2)
    deflections = mc_results['deflections'] * 1000  # 转换为mm
    ax2.hist(deflections, bins=30, density=True, alpha=0.7, color='coral', edgecolor='black')
    
    mu, std = norm.fit(deflections)
    x = np.linspace(deflections.min(), deflections.max(), 100)
    ax2.plot(x, norm.pdf(x, mu, std), 'r-', linewidth=2, label=f'正态拟合\nμ={mu:.3f}, σ={std:.3f}')
    
    ax2.axvline(mu, color='g', linestyle='--', linewidth=2, label=f'均值={mu:.3f}')
    ax2.axvline(mu - 2*std, color='orange', linestyle=':', linewidth=1.5, alpha=0.7)
    ax2.axvline(mu + 2*std, color='orange', linestyle=':', linewidth=1.5, alpha=0.7, label='95%置信区间')
    
    ax2.set_xlabel('中点挠度 (mm)', fontsize=11)
    ax2.set_ylabel('概率密度', fontsize=11)
    ax2.set_title('挠度概率分布', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 散点图矩阵
    ax3 = fig.add_subplot(1, 3, 3)
    
    E_samples = np.array(mc_results['param_samples']['E']) / 1e9
    I_samples = np.array(mc_results['param_samples']['I'])
    
    scatter = ax3.scatter(E_samples, I_samples, c=freqs, cmap='viridis', alpha=0.5, s=20)
    ax3.set_xlabel('弹性模量 E (GPa)', fontsize=11)
    ax3.set_ylabel('惯性矩 I (m^4)', fontsize=11)
    ax3.set_title('参数空间与响应关系', fontsize=12, fontweight='bold')
    cbar = plt.colorbar(scatter, ax=ax3)
    cbar.set_label('基频 (Hz)', fontsize=10)
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('案例2_不确定性量化.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("\n不确定性量化图已保存到: 案例2_不确定性量化.png")
    plt.close()


def visualize_sensitivity(sensitivities):
    """可视化敏感性分析结果"""
    fig = plt.figure(figsize=(16, 5))
    
    param_names = ['E', 'I', 'rho']
    param_labels = ['弹性模量 E (GPa)', '惯性矩 I (m^4)', '密度 ρ (kg/m³)']
    
    # 子图1: 参数对频率的影响
    ax1 = fig.add_subplot(1, 3, 1)
    colors = ['steelblue', 'coral', 'forestgreen']
    
    for i, (param, label, color) in enumerate(zip(param_names, param_labels, colors)):
        data = sensitivities[param]
        values = data['values']
        if param == 'E':
            values = values / 1e9
        elif param == 'rho':
            values = values
        
        ax1.plot(values, data['frequencies'], linewidth=2, label=label, color=color)
    
    ax1.set_xlabel('参数值', fontsize=11)
    ax1.set_ylabel('基频 (Hz)', fontsize=11)
    ax1.set_title('参数对频率的影响', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 参数对挠度的影响
    ax2 = fig.add_subplot(1, 3, 2)
    
    for i, (param, label, color) in enumerate(zip(param_names, param_labels, colors)):
        data = sensitivities[param]
        values = data['values']
        if param == 'E':
            values = values / 1e9
        
        ax2.plot(values, data['deflections']*1000, linewidth=2, label=label, color=color)
    
    ax2.set_xlabel('参数值', fontsize=11)
    ax2.set_ylabel('中点挠度 (mm)', fontsize=11)
    ax2.set_title('参数对挠度的影响', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 敏感性对比
    ax3 = fig.add_subplot(1, 3, 3)
    
    freq_sens = [sensitivities[p]['freq_sensitivity']*100 for p in param_names]
    deflection_sens = [sensitivities[p]['deflection_sensitivity']*100 for p in param_names]
    
    x = np.arange(len(param_names))
    width = 0.35
    
    bars1 = ax3.bar(x - width/2, freq_sens, width, label='频率敏感性', color='steelblue', edgecolor='black')
    bars2 = ax3.bar(x + width/2, deflection_sens, width, label='挠度敏感性', color='coral', edgecolor='black')
    
    ax3.set_ylabel('敏感性系数 (%)', fontsize=11)
    ax3.set_title('参数敏感性对比', fontsize=12, fontweight='bold')
    ax3.set_xticks(x)
    ax3.set_xticklabels(param_names)
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bars in [bars1, bars2]:
        for bar in bars:
            height = bar.get_height()
            ax3.text(bar.get_x() + bar.get_width()/2., height,
                    f'{height:.1f}%', ha='center', va='bottom', fontsize=9)
    
    plt.tight_layout()
    plt.savefig('案例2_敏感性分析.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("敏感性分析图已保存到: 案例2_敏感性分析.png")
    plt.close()


def visualize_validation(mc_results, freq_exp, deflection_exp):
    """可视化模型确认结果"""
    fig = plt.figure(figsize=(16, 5))
    
    # 子图1: 频率对比
    ax1 = fig.add_subplot(1, 3, 1)
    
    freq_sim = mc_results['frequencies']
    
    # 箱线图
    bp = ax1.boxplot([freq_sim, freq_exp], labels=['仿真', '试验'], patch_artist=True)
    bp['boxes'][0].set_facecolor('steelblue')
    bp['boxes'][1].set_facecolor('coral')
    
    ax1.set_ylabel('基频 (Hz)', fontsize=11)
    ax1.set_title('频率: 仿真 vs 试验', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3, axis='y')
    
    # 添加统计信息
    ax1.text(0.02, 0.98, f'仿真: μ={np.mean(freq_sim):.3f}, σ={np.std(freq_sim):.3f}\n'
                         f'试验: μ={np.mean(freq_exp):.3f}, σ={np.std(freq_exp):.3f}',
            transform=ax1.transAxes, fontsize=9, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    
    # 子图2: 挠度对比
    ax2 = fig.add_subplot(1, 3, 2)
    
    deflection_sim = mc_results['deflections'] * 1000
    deflection_exp_mm = deflection_exp * 1000
    
    bp = ax2.boxplot([deflection_sim, deflection_exp_mm], labels=['仿真', '试验'], patch_artist=True)
    bp['boxes'][0].set_facecolor('steelblue')
    bp['boxes'][1].set_facecolor('coral')
    
    ax2.set_ylabel('中点挠度 (mm)', fontsize=11)
    ax2.set_title('挠度: 仿真 vs 试验', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3, axis='y')
    
    ax2.text(0.02, 0.98, f'仿真: μ={np.mean(deflection_sim):.3f}, σ={np.std(deflection_sim):.3f}\n'
                         f'试验: μ={np.mean(deflection_exp_mm):.3f}, σ={np.std(deflection_exp_mm):.3f}',
            transform=ax2.transAxes, fontsize=9, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    
    # 子图3: 相关性分析
    ax3 = fig.add_subplot(1, 3, 3)
    
    # 使用仿真均值作为预测值
    freq_sim_mean = np.mean(freq_sim)
    deflection_sim_mean = np.mean(deflection_sim)
    
    # 散点图
    ax3.scatter(freq_exp, [freq_sim_mean]*len(freq_exp), 
               alpha=0.6, s=100, c='steelblue', label='频率', edgecolors='black')
    ax3.scatter(deflection_exp_mm, [deflection_sim_mean]*len(deflection_exp_mm), 
               alpha=0.6, s=100, c='coral', label='挠度', edgecolors='black', marker='s')
    
    # 完美预测线
    all_data = np.concatenate([freq_exp, deflection_exp_mm])
    min_val, max_val = all_data.min(), all_data.max()
    ax3.plot([min_val, max_val], [min_val, max_val], 'k--', linewidth=2, label='完美预测', alpha=0.5)
    
    ax3.set_xlabel('试验值', fontsize=11)
    ax3.set_ylabel('仿真值', fontsize=11)
    ax3.set_title('仿真-试验相关性', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('案例2_模型确认.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("模型确认图已保存到: 案例2_模型确认.png")
    plt.close()


def create_uncertainty_animation(mc_results):
    """创建不确定性传播动画"""
    print("\n正在生成不确定性动画...")
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    freqs = mc_results['frequencies']
    deflections = mc_results['deflections'] * 1000
    
    # 左图: 频率分布演化
    ax1 = axes[0]
    ax1.set_xlim(freqs.min()*0.9, freqs.max()*1.1)
    ax1.set_ylim(0, 50)
    ax1.set_xlabel('基频 (Hz)', fontsize=11)
    ax1.set_ylabel('频数', fontsize=11)
    ax1.set_title('基频分布收敛过程', fontsize=13, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 右图: 挠度分布演化
    ax2 = axes[1]
    ax2.set_xlim(deflections.min()*0.9, deflections.max()*1.1)
    ax2.set_ylim(0, 50)
    ax2.set_xlabel('中点挠度 (mm)', fontsize=11)
    ax2.set_ylabel('频数', fontsize=11)
    ax2.set_title('挠度分布收敛过程', fontsize=13, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    def init():
        return []
    
    def update(frame):
        ax1.clear()
        ax2.clear()
        
        n_samples = (frame + 1) * 50
        if n_samples > len(freqs):
            n_samples = len(freqs)
        
        # 重新设置轴属性
        ax1.set_xlim(freqs.min()*0.9, freqs.max()*1.1)
        ax1.set_ylim(0, 50)
        ax1.set_xlabel('基频 (Hz)', fontsize=11)
        ax1.set_ylabel('频数', fontsize=11)
        ax1.set_title(f'基频分布 (n={n_samples})', fontsize=13, fontweight='bold')
        ax1.grid(True, alpha=0.3)
        
        ax2.set_xlim(deflections.min()*0.9, deflections.max()*1.1)
        ax2.set_ylim(0, 50)
        ax2.set_xlabel('中点挠度 (mm)', fontsize=11)
        ax2.set_ylabel('频数', fontsize=11)
        ax2.set_title(f'挠度分布 (n={n_samples})', fontsize=13, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        
        # 绘制直方图
        ax1.hist(freqs[:n_samples], bins=20, alpha=0.7, color='steelblue', edgecolor='black')
        ax2.hist(deflections[:n_samples], bins=20, alpha=0.7, color='coral', edgecolor='black')
        
        # 添加均值线
        ax1.axvline(np.mean(freqs[:n_samples]), color='r', linestyle='--', linewidth=2)
        ax2.axvline(np.mean(deflections[:n_samples]), color='r', linestyle='--', linewidth=2)
        
        return []
    
    n_frames = len(freqs) // 50
    anim = FuncAnimation(fig, update, init_func=init, frames=n_frames,
                        interval=100, blit=False)
    
    writer = PillowWriter(fps=10)
    anim.save('案例2_不确定性收敛.gif', writer=writer)
    print("动画已保存到: 案例2_不确定性收敛.gif")
    plt.close()


if __name__ == "__main__":
    main()

案例3: 仿真结果与试验数据对比验证

目标：确认仿真模型对真实结构的预测能力

方法：

1. 模态试验对比
2. 静载试验对比
3. 动态响应对比

预期成果：

• 计算确认度量
• 评估模型适用性
• 建立置信区间

import matplotlib
matplotlib.use('Agg')
"""
主题097 - 案例3: 仿真结果与试验数据对比验证
对比分析仿真预测与物理试验结果，评估模型可信度
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.patches import FancyBboxPatch, Rectangle, Circle, FancyArrowPatch
from scipy import stats
from scipy.stats import norm, ttest_ind, pearsonr

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


def generate_simulation_data(n_points=100):
    """生成仿真数据"""
    np.random.seed(42)
    
    # 频率响应函数 (FRF)
    freq = np.linspace(0, 50, n_points)  # 0-50 Hz
    
    # 理论模态频率
    f1, f2, f3 = 5.0, 15.0, 30.0
    
    # 仿真FRF (包含三个模态)
    H_sim = (1 / (1 - (freq/f1)**2 + 0.05j*(freq/f1))) + \
            (0.8 / (1 - (freq/f2)**2 + 0.03j*(freq/f2))) + \
            (0.6 / (1 - (freq/f3)**2 + 0.04j*(freq/f3)))
    
    # 添加数值噪声
    noise = 0.02 * (np.random.randn(n_points) + 1j * np.random.randn(n_points))
    H_sim += noise
    
    return freq, H_sim


def generate_experimental_data(freq, H_sim, noise_level=0.1):
    """生成试验数据（在仿真数据基础上添加试验噪声）"""
    np.random.seed(123)
    
    # 添加试验噪声（比仿真噪声大）
    noise_real = noise_level * np.random.randn(len(freq))
    noise_imag = noise_level * np.random.randn(len(freq))
    
    H_exp = H_sim + (noise_real + 1j * noise_imag)
    
    # 添加系统误差（频率偏移）
    freq_shift = 0.02  # 2%频率偏移
    H_exp_shifted = np.interp(freq * (1 + freq_shift), freq, H_exp.real) + \
                   1j * np.interp(freq * (1 + freq_shift), freq, H_exp.imag)
    
    return H_exp_shifted


def calculate_frf_metrics(freq, H_sim, H_exp):
    """计算FRF对比指标"""
    
    # 幅值误差
    mag_sim = np.abs(H_sim)
    mag_exp = np.abs(H_exp)
    
    # 相位误差
    phase_sim = np.angle(H_sim, deg=True)
    phase_exp = np.angle(H_exp, deg=True)
    
    # 幅值均方根误差
    rmse_mag = np.sqrt(np.mean((mag_sim - mag_exp)**2))
    
    # 相位均方根误差
    phase_diff = np.unwrap(np.angle(H_sim)) - np.unwrap(np.angle(H_exp))
    rmse_phase = np.sqrt(np.mean(phase_diff**2)) * 180 / np.pi
    
    # 相关系数
    corr_mag, p_mag = pearsonr(mag_sim, mag_exp)
    
    # 频率响应误差 (复数)
    frf_error = np.abs(H_sim - H_exp)
    mean_frf_error = np.mean(frf_error)
    max_frf_error = np.max(frf_error)
    
    # 模态参数提取（简化）
    peaks_sim = detect_peaks(mag_sim)
    peaks_exp = detect_peaks(mag_exp)
    
    return {
        'RMSE_Magnitude': rmse_mag,
        'RMSE_Phase': rmse_phase,
        'Correlation': corr_mag,
        'P_Value': p_mag,
        'Mean_FRF_Error': mean_frf_error,
        'Max_FRF_Error': max_frf_error,
        'Peaks_Sim': peaks_sim,
        'Peaks_Exp': peaks_exp
    }


def detect_peaks(magnitude, threshold=0.5):
    """简单的峰值检测"""
    peaks = []
    for i in range(1, len(magnitude) - 1):
        if magnitude[i] > magnitude[i-1] and magnitude[i] > magnitude[i+1]:
            if magnitude[i] > threshold * np.max(magnitude):
                peaks.append(i)
    return peaks


def time_domain_comparison():
    """时域响应对比分析"""
    print("\n" + "="*60)
    print("时域响应对比分析")
    print("="*60)
    
    # 时间向量
    t = np.linspace(0, 10, 1000)
    dt = t[1] - t[0]
    
    # 激励信号 (扫频)
    f0, f1 = 1, 20
    excitation = np.sin(2 * np.pi * (f0 + (f1 - f0) * t / 10) * t)
    
    # 系统参数
    m, c, k = 1.0, 0.5, 100.0
    omega_n = np.sqrt(k/m)
    zeta = c / (2 * np.sqrt(m*k))
    
    # 仿真响应 (精确解)
    response_sim = simulate_sdof_response(m, c, k, excitation, dt)
    
    # 试验响应 (添加噪声和模型误差)
    np.random.seed(456)
    
    # 模型误差：质量偏差5%
    m_err = m * 1.05
    response_exp = simulate_sdof_response(m_err, c, k, excitation, dt)
    
    # 添加测量噪声
    noise = 0.1 * np.max(np.abs(response_exp)) * np.random.randn(len(t))
    response_exp += noise
    
    # 计算误差指标
    error = response_sim - response_exp
    rmse = np.sqrt(np.mean(error**2))
    mae = np.mean(np.abs(error))
    max_error = np.max(np.abs(error))
    
    print(f"\n时域响应误差指标:")
    print(f"  RMSE: {rmse:.6f}")
    print(f"  MAE: {mae:.6f}")
    print(f"  最大误差: {max_error:.6f}")
    
    # 计算相关系数
    corr, p_value = pearsonr(response_sim, response_exp)
    print(f"  相关系数: {corr:.4f} (p={p_value:.4e})")
    
    return {
        't': t,
        'excitation': excitation,
        'response_sim': response_sim,
        'response_exp': response_exp,
        'error': error,
        'metrics': {
            'RMSE': rmse,
            'MAE': mae,
            'MaxError': max_error,
            'Correlation': corr
        }
    }


def simulate_sdof_response(m, c, k, force, dt):
    """模拟单自由度系统响应（Newmark-beta方法）"""
    n = len(force)
    u = np.zeros(n)
    v = np.zeros(n)
    a = np.zeros(n)
    
    # Newmark参数
    beta = 0.25
    gamma = 0.5
    
    # 初始条件
    a[0] = (force[0] - c*v[0] - k*u[0]) / m
    
    for i in range(n - 1):
        # 预测
        u_pred = u[i] + dt*v[i] + (dt**2/2)*(1-2*beta)*a[i]
        v_pred = v[i] + dt*(1-gamma)*a[i]
        
        # 求解加速度
        a[i+1] = (force[i+1] - c*v_pred - k*u_pred) / (m + c*dt*gamma + k*dt**2*beta)
        
        # 校正
        u[i+1] = u_pred + dt**2*beta*a[i+1]
        v[i+1] = v_pred + dt*gamma*a[i+1]
    
    return u


def statistical_validation():
    """统计验证分析"""
    print("\n" + "="*60)
    print("统计验证分析")
    print("="*60)
    
    # 生成多组试验数据
    np.random.seed(789)
    n_samples = 50
    
    # 仿真结果（理论值）
    sim_results = {
        'frequency': 5.0,
        'damping': 0.02,
        'amplitude': 1.0
    }
    
    # 试验结果（添加随机误差）
    exp_frequencies = sim_results['frequency'] * (1 + 0.03 * np.random.randn(n_samples))
    exp_dampings = sim_results['damping'] * (1 + 0.15 * np.random.randn(n_samples))
    exp_amplitudes = sim_results['amplitude'] * (1 + 0.05 * np.random.randn(n_samples))
    
    # 统计检验
    print("\n模态频率验证:")
    print(f"  仿真值: {sim_results['frequency']:.3f} Hz")
    print(f"  试验均值: {np.mean(exp_frequencies):.3f} Hz")
    print(f"  试验标准差: {np.std(exp_frequencies):.3f} Hz")
    
    # t检验（单样本t检验）
    from scipy.stats import ttest_1samp
    t_stat, p_val = ttest_1samp(exp_frequencies, sim_results['frequency'])
    print(f"  t统计量: {t_stat:.3f}, p值: {p_val:.4f}")
    
    if p_val > 0.05:
        print("  结论: 仿真与试验无显著差异 (p > 0.05)")
    else:
        print("  结论: 仿真与试验存在显著差异 (p < 0.05)")
    
    print("\n阻尼比验证:")
    print(f"  仿真值: {sim_results['damping']:.4f}")
    print(f"  试验均值: {np.mean(exp_dampings):.4f}")
    print(f"  试验标准差: {np.std(exp_dampings):.4f}")
    
    print("\n振幅验证:")
    print(f"  仿真值: {sim_results['amplitude']:.3f}")
    print(f"  试验均值: {np.mean(exp_amplitudes):.3f}")
    print(f"  试验标准差: {np.std(exp_amplitudes):.3f}")
    
    return {
        'sim': sim_results,
        'exp_freq': exp_frequencies,
        'exp_damp': exp_dampings,
        'exp_amp': exp_amplitudes
    }


def main():
    """主函数"""
    print("=" * 70)
    print("主题097 - 案例3: 仿真结果与试验数据对比验证")
    print("=" * 70)
    
    # 1. 频域对比分析
    print("\n" + "="*60)
    print("1. 频域响应函数(FRF)对比")
    print("="*60)
    
    freq, H_sim = generate_simulation_data(n_points=200)
    H_exp = generate_experimental_data(freq, H_sim, noise_level=0.15)
    
    frf_metrics = calculate_frf_metrics(freq, H_sim, H_exp)
    
    print("\nFRF对比指标:")
    for key, value in frf_metrics.items():
        if not isinstance(value, list):
            print(f"  {key}: {value:.4f}")
    
    # 2. 时域对比分析
    time_data = time_domain_comparison()
    
    # 3. 统计验证
    stat_data = statistical_validation()
    
    # 生成可视化
    print("\n" + "="*60)
    print("生成可视化结果")
    print("="*60)
    
    visualize_frf_comparison(freq, H_sim, H_exp, frf_metrics)
    visualize_time_comparison(time_data)
    visualize_statistical_validation(stat_data)
    create_validation_animation(freq, H_sim, H_exp)
    
    print("\n" + "=" * 70)
    print("案例3完成！所有结果已保存。")
    print("=" * 70)


def visualize_frf_comparison(freq, H_sim, H_exp, metrics):
    """可视化FRF对比结果"""
    fig = plt.figure(figsize=(16, 10))
    
    mag_sim = np.abs(H_sim)
    mag_exp = np.abs(H_exp)
    phase_sim = np.angle(H_sim, deg=True)
    phase_exp = np.angle(H_exp, deg=True)
    
    # 子图1: 幅值对比
    ax1 = fig.add_subplot(2, 2, 1)
    ax1.semilogy(freq, mag_sim, 'b-', linewidth=2, label='仿真', alpha=0.8)
    ax1.semilogy(freq, mag_exp, 'r--', linewidth=2, label='试验', alpha=0.8)
    
    # 标记峰值
    for peak in metrics['Peaks_Sim'][:3]:
        ax1.plot(freq[peak], mag_sim[peak], 'bo', markersize=10)
    for peak in metrics['Peaks_Exp'][:3]:
        ax1.plot(freq[peak], mag_exp[peak], 'rs', markersize=8)
    
    ax1.set_xlabel('频率 (Hz)', fontsize=11)
    ax1.set_ylabel('幅值', fontsize=11)
    ax1.set_title('FRF幅值对比', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3, which='both')
    ax1.set_xlim(0, 50)
    
    # 子图2: 相位对比
    ax2 = fig.add_subplot(2, 2, 2)
    ax2.plot(freq, phase_sim, 'b-', linewidth=2, label='仿真', alpha=0.8)
    ax2.plot(freq, phase_exp, 'r--', linewidth=2, label='试验', alpha=0.8)
    ax2.set_xlabel('频率 (Hz)', fontsize=11)
    ax2.set_ylabel('相位 (度)', fontsize=11)
    ax2.set_title('FRF相位对比', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 50)
    
    # 子图3: 幅值误差
    ax3 = fig.add_subplot(2, 2, 3)
    mag_error = np.abs(mag_sim - mag_exp)
    ax3.fill_between(freq, mag_error, alpha=0.5, color='orange')
    ax3.plot(freq, mag_error, 'g-', linewidth=1.5)
    ax3.axhline(y=metrics['RMSE_Magnitude'], color='r', linestyle='--', 
               linewidth=2, label=f'RMSE = {metrics["RMSE_Magnitude"]:.3f}')
    ax3.set_xlabel('频率 (Hz)', fontsize=11)
    ax3.set_ylabel('幅值误差', fontsize=11)
    ax3.set_title('幅值误差分布', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=10)
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim(0, 50)
    
    # 子图4: 复平面表示
    ax4 = fig.add_subplot(2, 2, 4)
    ax4.plot(H_sim.real, H_sim.imag, 'b-', linewidth=2, label='仿真', alpha=0.7)
    ax4.plot(H_exp.real, H_exp.imag, 'r--', linewidth=2, label='试验', alpha=0.7)
    ax4.set_xlabel('实部', fontsize=11)
    ax4.set_ylabel('虚部', fontsize=11)
    ax4.set_title('FRF复平面表示', fontsize=12, fontweight='bold')
    ax4.legend(fontsize=10)
    ax4.grid(True, alpha=0.3)
    ax4.axis('equal')
    
    plt.tight_layout()
    plt.savefig('案例3_FRF对比分析.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("\nFRF对比分析图已保存到: 案例3_FRF对比分析.png")
    plt.close()


def visualize_time_comparison(time_data):
    """可视化时域对比结果"""
    fig = plt.figure(figsize=(16, 5))
    
    t = time_data['t']
    excitation = time_data['excitation']
    response_sim = time_data['response_sim']
    response_exp = time_data['response_exp']
    error = time_data['error']
    
    # 子图1: 激励信号
    ax1 = fig.add_subplot(1, 3, 1)
    ax1.plot(t, excitation, 'k-', linewidth=1.5)
    ax1.set_xlabel('时间 (s)', fontsize=11)
    ax1.set_ylabel('激励幅值', fontsize=11)
    ax1.set_title('激励信号（扫频）', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, 10)
    
    # 子图2: 响应对比
    ax2 = fig.add_subplot(1, 3, 2)
    ax2.plot(t, response_sim, 'b-', linewidth=2, label='仿真', alpha=0.8)
    ax2.plot(t, response_exp, 'r--', linewidth=2, label='试验', alpha=0.8)
    ax2.set_xlabel('时间 (s)', fontsize=11)
    ax2.set_ylabel('位移响应', fontsize=11)
    ax2.set_title('时域响应对比', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 10)
    
    # 子图3: 误差分析
    ax3 = fig.add_subplot(1, 3, 3)
    ax3.fill_between(t, error, alpha=0.5, color='orange', label='误差')
    ax3.plot(t, error, 'g-', linewidth=1.5)
    ax3.axhline(y=0, color='k', linestyle='-', linewidth=1)
    
    # 添加误差指标
    metrics = time_data['metrics']
    ax3.axhline(y=metrics['RMSE'], color='r', linestyle='--', 
               linewidth=2, label=f'RMSE = {metrics["RMSE"]:.4f}')
    ax3.axhline(y=-metrics['RMSE'], color='r', linestyle='--', linewidth=2)
    
    ax3.set_xlabel('时间 (s)', fontsize=11)
    ax3.set_ylabel('误差', fontsize=11)
    ax3.set_title('响应误差', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim(0, 10)
    
    plt.tight_layout()
    plt.savefig('案例3_时域对比分析.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("时域对比分析图已保存到: 案例3_时域对比分析.png")
    plt.close()


def visualize_statistical_validation(stat_data):
    """可视化统计验证结果"""
    fig = plt.figure(figsize=(16, 5))
    
    # 子图1: 频率分布
    ax1 = fig.add_subplot(1, 3, 1)
    
    exp_freq = stat_data['exp_freq']
    sim_freq = stat_data['sim']['frequency']
    
    ax1.hist(exp_freq, bins=15, density=True, alpha=0.7, color='steelblue', 
            edgecolor='black', label='试验分布')
    ax1.axvline(sim_freq, color='r', linestyle='--', linewidth=3, 
               label=f'仿真值 = {sim_freq:.2f} Hz')
    ax1.axvline(np.mean(exp_freq), color='g', linestyle='-', linewidth=2, 
               label=f'试验均值 = {np.mean(exp_freq):.2f} Hz')
    
    # 添加置信区间
    ci_low = np.percentile(exp_freq, 2.5)
    ci_high = np.percentile(exp_freq, 97.5)
    ax1.axvspan(ci_low, ci_high, alpha=0.2, color='orange', label='95% CI')
    
    ax1.set_xlabel('频率 (Hz)', fontsize=11)
    ax1.set_ylabel('概率密度', fontsize=11)
    ax1.set_title('模态频率统计分布', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 阻尼比分布
    ax2 = fig.add_subplot(1, 3, 2)
    
    exp_damp = stat_data['exp_damp']
    sim_damp = stat_data['sim']['damping']
    
    ax2.hist(exp_damp, bins=15, density=True, alpha=0.7, color='coral', 
            edgecolor='black', label='试验分布')
    ax2.axvline(sim_damp, color='r', linestyle='--', linewidth=3, 
               label=f'仿真值 = {sim_damp:.3f}')
    ax2.axvline(np.mean(exp_damp), color='g', linestyle='-', linewidth=2, 
               label=f'试验均值 = {np.mean(exp_damp):.3f}')
    
    ax2.set_xlabel('阻尼比', fontsize=11)
    ax2.set_ylabel('概率密度', fontsize=11)
    ax2.set_title('阻尼比统计分布', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 振幅分布
    ax3 = fig.add_subplot(1, 3, 3)
    
    exp_amp = stat_data['exp_amp']
    sim_amp = stat_data['sim']['amplitude']
    
    ax3.hist(exp_amp, bins=15, density=True, alpha=0.7, color='forestgreen', 
            edgecolor='black', label='试验分布')
    ax3.axvline(sim_amp, color='r', linestyle='--', linewidth=3, 
               label=f'仿真值 = {sim_amp:.2f}')
    ax3.axvline(np.mean(exp_amp), color='g', linestyle='-', linewidth=2, 
               label=f'试验均值 = {np.mean(exp_amp):.2f}')
    
    ax3.set_xlabel('振幅', fontsize=11)
    ax3.set_ylabel('概率密度', fontsize=11)
    ax3.set_title('振幅统计分布', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('案例3_统计验证分析.png', dpi=150, bbox_inches='tight',
               facecolor='white', edgecolor='none')
    print("统计验证分析图已保存到: 案例3_统计验证分析.png")
    plt.close()


def create_validation_animation(freq, H_sim, H_exp):
    """创建验证过程动画"""
    print("\n正在生成验证动画...")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    mag_sim = np.abs(H_sim)
    mag_exp = np.abs(H_exp)
    
    # 设置轴属性
    for ax in axes.flat:
        ax.set_xlim(0, 50)
        ax.grid(True, alpha=0.3)
    
    axes[0, 0].set_ylabel('幅值', fontsize=11)
    axes[0, 0].set_title('FRF幅值对比', fontsize=12, fontweight='bold')
    axes[0, 0].set_yscale('log')
    
    axes[0, 1].set_ylabel('相位 (度)', fontsize=11)
    axes[0, 1].set_title('FRF相位对比', fontsize=12, fontweight='bold')
    
    axes[1, 0].set_xlabel('频率 (Hz)', fontsize=11)
    axes[1, 0].set_ylabel('幅值误差', fontsize=11)
    axes[1, 0].set_title('误差演化', fontsize=12, fontweight='bold')
    
    axes[1, 1].set_xlabel('频率 (Hz)', fontsize=11)
    axes[1, 1].set_ylabel('相关系数', fontsize=11)
    axes[1, 1].set_title('相关系数演化', fontsize=12, fontweight='bold')
    
    def init():
        return []
    
    def update(frame):
        # 清除之前的线条
        for ax in axes.flat:
            ax.clear()
            ax.set_xlim(0, 50)
            ax.grid(True, alpha=0.3)
        
        idx = (frame + 1) * 5
        if idx > len(freq):
            idx = len(freq)
        
        # 子图1: 幅值
        axes[0, 0].semilogy(freq[:idx], mag_sim[:idx], 'b-', linewidth=2, label='仿真')
        axes[0, 0].semilogy(freq[:idx], mag_exp[:idx], 'r--', linewidth=2, label='试验')
        axes[0, 0].set_ylabel('幅值', fontsize=11)
        axes[0, 0].set_title(f'FRF幅值对比 (f < {freq[idx-1]:.1f} Hz)', fontsize=12, fontweight='bold')
        axes[0, 0].legend(fontsize=9)
        
        # 子图2: 相位
        phase_sim = np.angle(H_sim[:idx], deg=True)
        phase_exp = np.angle(H_exp[:idx], deg=True)
        axes[0, 1].plot(freq[:idx], phase_sim, 'b-', linewidth=2, label='仿真')
        axes[0, 1].plot(freq[:idx], phase_exp, 'r--', linewidth=2, label='试验')
        axes[0, 1].set_ylabel('相位 (度)', fontsize=11)
        axes[0, 1].set_title('FRF相位对比', fontsize=12, fontweight='bold')
        axes[0, 1].legend(fontsize=9)
        
        # 子图3: 误差
        error = np.abs(mag_sim[:idx] - mag_exp[:idx])
        axes[1, 0].fill_between(freq[:idx], error, alpha=0.5, color='orange')
        axes[1, 0].plot(freq[:idx], error, 'g-', linewidth=1.5)
        axes[1, 0].set_xlabel('频率 (Hz)', fontsize=11)
        axes[1, 0].set_ylabel('幅值误差', fontsize=11)
        axes[1, 0].set_title('误差演化', fontsize=12, fontweight='bold')
        
        # 子图4: 累积相关系数
        if idx > 10:
            corr = pearsonr(mag_sim[:idx], mag_exp[:idx])[0]
        else:
            corr = 0
        
        axes[1, 1].bar([1], [corr], color='steelblue', alpha=0.7, width=0.5)
        axes[1, 1].set_ylim(0, 1)
        axes[1, 1].set_xlabel('验证指标', fontsize=11)
        axes[1, 1].set_ylabel('相关系数', fontsize=11)
        axes[1, 1].set_title(f'相关系数: {corr:.3f}', fontsize=12, fontweight='bold')
        axes[1, 1].set_xticks([1])
        axes[1, 1].set_xticklabels(['幅值相关'])
        
        return []
    
    n_frames = len(freq) // 5
    anim = FuncAnimation(fig, update, init_func=init, frames=n_frames,
                        interval=100, blit=False)
    
    writer = PillowWriter(fps=10)
    anim.save('案例3_验证过程动画.gif', writer=writer)
    print("动画已保存到: 案例3_验证过程动画.gif")
    plt.close()


if __name__ == "__main__":
    main()

---

总结与展望

V&V最佳实践

1. 早期规划：在项目启动阶段制定V&V计划
2. 文档化：详细记录所有验证和确认活动
3. 层次化方法：从简单到复杂逐步验证
4. 不确定性量化：始终伴随不确定度估计
5. 持续更新：根据新数据更新模型和确认状态

未来发展趋势

1. 自动化V&V：基于机器学习的自动验证工具
2. 实时确认：数字孪生中的在线模型确认
3. 多尺度V&V：跨尺度模型的统一验证框架
4. 云仿真验证：分布式仿真环境的验证方法

---