主题082：深度学习求解辐射传递方程

摘要

辐射传递方程（RTE）的求解是辐射换热计算的核心问题，传统数值方法如离散坐标法、有限体积法等计算成本高，难以满足实时预测需求。物理信息神经网络（PINN）将物理定律嵌入神经网络，为求解RTE提供了全新途径。本主题系统介绍如何利用深度学习方法直接求解辐射传递方程，包括PINN的基本原理、网络架构设计、损失函数构建和训练策略。通过将RTE作为约束嵌入神经网络，实现辐射场的快速准确预测，为实时辐射换热分析和优化设计提供强大工具。

关键词

物理信息神经网络、PINN、辐射传递方程、深度学习、实时求解、代理模型、物理约束、自动微分

---

1. 引言：辐射传递方程求解的挑战

1.1 辐射传递方程简介

辐射传递方程（Radiative Transfer Equation, RTE）描述了辐射能量在介质中的传播、吸收、发射和散射过程。对于一维灰体介质，RTE可以表示为：

$$\mu \frac{dI}{d\tau }+I=(1-\omega )I_b+\frac{\omega }{4\pi }{\int }_{4\pi }I\Phi d\Omega$$

其中：

• $I$ 是辐射强度
• $\tau$ 是光学厚度
• $\mu$ 是方向余弦
• $\omega$ 是散射反照率
• $I_b$ 是黑体辐射强度
• $\Phi$ 是散射相函数

1.2 传统数值方法的局限性

计算成本高

传统方法如离散坐标法（DOM）、有限体积法（FVM）需要离散空间和方向：

• 空间离散：$N_x\times N_y\times N_z$ 个网格
• 方向离散：$N_{\theta }\times N_{\varphi }$ 个方向
• 总自由度：$N_x\times N_y\times N_z\times N_{\theta }\times N_{\varphi }$

对于复杂三维问题，计算量巨大。

实时性差

• 每个工况需要重新求解
• 优化设计需要成千上万次求解
• 难以满足实时控制需求

网格依赖

• 需要精细网格保证精度
• 网格生成耗时
• 复杂几何处理困难

1.3 深度学习的优势

快速推理

训练后的神经网络预测速度极快（毫秒级），适合实时应用。

无网格方法

神经网络是连续函数近似，不需要空间离散。

物理一致性

PINN将物理方程作为约束，保证预测结果满足物理定律。

泛化能力

学习到的解可以推广到新的参数和边界条件。

---

2. 物理信息神经网络（PINN）基础

2.1 PINN基本原理

物理信息神经网络（Physics-Informed Neural Networks, PINN）由Raissi等人于2019年提出，核心思想是将物理定律嵌入神经网络的损失函数。

网络架构

PINN使用标准的全连接神经网络：

$$\hat{u}(\mathbf{x};\theta )=\mathcal{N}(\mathbf{x};\theta )$$

其中$\mathbf{x}$是输入（空间坐标、时间、参数等），$\theta$是网络参数。

损失函数

PINN的损失函数包含两部分：

$$\mathcal{L}={\mathcal{L}}_{data}+\lambda {\mathcal{L}}_{physics}$$

• 数据损失 ${\mathcal{L}}_{data}$：拟合观测数据
• 物理损失 ${\mathcal{L}}_{physics}$：满足物理方程
• 权重 $\lambda$：平衡两项的重要性

2.2 自动微分

PINN的关键技术是自动微分（Automatic Differentiation, AD），用于计算神经网络输出的导数。

前向模式

计算函数值的同时计算导数。

反向模式

先计算函数值，再反向传播计算梯度。这是深度学习框架（如PyTorch、TensorFlow）的标准方法。

应用于PDE

对于偏微分方程：

$$\mathcal{F}(u,u_x,u_{xx},...)=0$$

使用自动微分计算：

$$u_x=\frac{\partial \hat{u}}{\partial x},u_{xx}=\frac{{\partial }^2\hat{u}}{\partial x^2},...$$

2.3 物理损失构建

残差计算

将神经网络解代入PDE，计算残差：

$$r=\mathcal{F}(\hat{u},{\hat{u}}_x,{\hat{u}}_{xx},...)$$

损失函数

$${\mathcal{L}}_{physics}=\frac{1}{N}\sum _{i=1}^N|r({\mathbf{x}}_i){|}^2$$

边界条件

边界条件也作为约束：

$${\mathcal{L}}_{BC}=\frac{1}{N_{BC}}\sum _{i=1}^{N_{BC}}|\hat{u}({\mathbf{x}}_i^{BC})-u^{BC}({\mathbf{x}}_i^{BC}){|}^2$$

---

3. PINN求解辐射传递方程

3.1 RTE的PINN形式

对于一维灰体介质的RTE：

$$\mu \frac{\partial I}{\partial \tau }+I=S$$

其中源项 $S=(1-\omega )I_b+\frac{\omega }{4\pi }\int I\Phi d\Omega$。

神经网络近似

构建神经网络 $\hat{I}(\tau ,\mu ;\theta )$ 近似辐射强度。

残差定义

$$r_{RTE}=\mu \frac{\partial \hat{I}}{\partial \tau }+\hat{I}-\hat{S}$$

边界残差

在边界$\tau =0$和$\tau ={\tau }_L$：

$$r_{BC,0}=\hat{I}(0,\mu >0)-I_{in,0}$$

$$r_{BC,L}=\hat{I}({\tau }_L,\mu <0)-I_{in,L}$$

3.2 网络架构设计

输入层

• 光学厚度$\tau \in [0,{\tau }_L]$
• 方向余弦$\mu \in [-1,1]$
• 可选：散射反照率$\omega$、光学厚度${\tau }_L$等参数

隐藏层

• 层数：4-8层
• 每层节点数：50-200
• 激活函数：tanh或sin（适合周期性/振荡解）

输出层

• 辐射强度$I(\tau ,\mu )$
• 可以使用归一化输出

3.3 损失函数构建

总损失

$$\mathcal{L}={\lambda }_{RTE}{\mathcal{L}}_{RTE}+{\lambda }_{BC}{\mathcal{L}}_{BC}+{\lambda }_{data}{\mathcal{L}}_{data}$$

各项定义

$${\mathcal{L}}_{RTE}=\frac{1}{N_{col}}\sum _{i=1}^{N_{col}}|r_{RTE}({\tau }_i,{\mu }_i){|}^2$$

$${\mathcal{L}}_{BC}=\frac{1}{N_{BC}}\sum _{i=1}^{N_{BC}}|r_{BC}({\tau }_i^{BC},{\mu }_i^{BC}){|}^2$$

$${\mathcal{L}}_{data}=\frac{1}{N_{data}}\sum _{i=1}^{N_{data}}|\hat{I}({\tau }_i^{data},{\mu }_i^{data})-I_i^{data}{|}^2$$

权重选择

• ${\lambda }_{RTE}=1$（主要约束）
• ${\lambda }_{BC}=10-100$（边界条件更重要）
• ${\lambda }_{data}$：有数据时使用

---

4. 训练策略与优化

4.1 采样策略

配点采样

在计算域内随机采样点计算物理损失：

• 均匀采样
• 拉丁超立方采样（LHS）
• 自适应采样（在残差大的区域加密）

边界采样

在边界上采样点施加边界条件。

数据采样

如果有观测数据，在数据点位置采样。

4.2 训练技巧

学习率调度

• 初始学习率：1e-3
• 学习率衰减：每1000步衰减0.9倍
• 早停：验证损失不再下降时停止

权重自适应

使用自适应权重平衡不同损失项：

$${\lambda }_i=\frac{|{\nabla }_{\theta }{\mathcal{L}}_i|}{{\sum }_j|{\nabla }_{\theta }{\mathcal{L}}_j|}$$

课程学习

从简单问题开始，逐步增加难度。

4.3 收敛诊断

残差监控

监控各项损失的收敛情况。

物理一致性检查

验证预测解是否满足物理约束（如能量守恒）。

与解析解对比

对于简单问题，与解析解对比验证。

---

5. 扩展与应用

5.1 多维RTE

对于二维/三维问题：

$$\mathbf{s}\cdot \nabla I+\beta I=\kappa I_b+\frac{{\sigma }_s}{4\pi }\int I\Phi d\Omega$$

输入包括空间坐标$(x,y,z)$和方向$(\theta ,\varphi )$。

5.2 非灰体介质

考虑光谱依赖性：

$$\mu \frac{\partial I_{\nu }}{\partial {\tau }_{\nu }}+I_{\nu }=S_{\nu }$$

输入增加波数$\nu$或波长$\lambda$。

5.3 瞬态问题

加入时间维度：

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mu \frac{\partial I}{\partial \tau }+I=S$$

输入增加时间$t$。

5.4 逆问题求解

PINN可以同时求解正问题和逆问题：

• 正问题：已知物性，求解辐射场
• 逆问题：已知辐射场测量，反演物性参数

损失函数同时包含正问题和逆问题的约束。

---

6. Python仿真案例

案例1：PINN求解一维RTE

使用PINN求解一维灰体介质的辐射传递方程，与解析解对比。

案例2：不同散射条件下的求解

测试各向同性、前向和后向散射对PINN求解的影响。

案例3：PINN与传统方法对比

比较PINN与离散坐标法的精度和计算效率。

案例4：参数化PINN

训练一个可以处理不同光学厚度和散射反照率的通用模型。

案例5：逆问题求解

利用PINN从辐射场测量反演介质吸收系数。

案例6：PINN求解RTE的动态演示

创建动画展示PINN训练过程和预测结果。

---

附录：Python代码说明

本主题配套的Python程序run_simulation.py实现了6个PINN求解辐射传递方程的仿真案例：

依赖库

• numpy：数值计算
• matplotlib：数据可视化
• scipy：科学计算（提供解析解对比）

运行方法

python run_simulation.py

输出文件

• case1_pinn_basic.png：PINN基本求解结果
• case2_scattering_effects.png：不同散射条件对比
• case3_pinn_vs_dom.png：PINN与传统方法对比
• case4_parametric_pinn.png：参数化PINN结果
• case5_inverse_problem.png：逆问题求解结果
• pinn_rte_animation.gif：PINN求解RTE动态演示

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
主题082：深度学习求解辐射传递方程 - Python仿真程序
包含6个案例：
1. PINN求解一维RTE
2. 不同散射条件下的求解
3. PINN与传统方法对比
4. 参数化PINN
5. 逆问题求解
6. PINN求解RTE的动态演示

注意：本程序使用简化的PINN实现，仅用于教学演示
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyArrowPatch
import matplotlib.animation as animation
from scipy.integrate import solve_bvp
import warnings
warnings.filterwarnings('ignore')

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


class SimplePINN:
    """
    简化的物理信息神经网络（PINN）
    使用numpy实现基本的前馈神经网络
    """
    def __init__(self, layers, activation='tanh'):
        """
        初始化网络
        layers: 各层节点数列表，如[2, 50, 50, 1]
        activation: 激活函数类型
        """
        self.layers = layers
        self.n_layers = len(layers)
        self.activation = activation
        
        # 初始化权重和偏置
        self.weights = []
        self.biases = []
        for i in range(self.n_layers - 1):
            W = np.random.randn(layers[i], layers[i+1]) * np.sqrt(2.0 / layers[i])
            b = np.zeros((1, layers[i+1]))
            self.weights.append(W)
            self.biases.append(b)
    
    def activate(self, z):
        """激活函数"""
        if self.activation == 'tanh':
            return np.tanh(z)
        elif self.activation == 'relu':
            return np.maximum(0, z)
        elif self.activation == 'sigmoid':
            return 1 / (1 + np.exp(-z))
        else:
            return np.tanh(z)
    
    def activate_derivative(self, z):
        """激活函数的导数"""
        if self.activation == 'tanh':
            return 1 - np.tanh(z)**2
        elif self.activation == 'relu':
            return (z > 0).astype(float)
        elif self.activation == 'sigmoid':
            s = 1 / (1 + np.exp(-z))
            return s * (1 - s)
        else:
            return 1 - np.tanh(z)**2
    
    def forward(self, x):
        """
        前向传播
        x: 输入，形状为 (n_samples, n_inputs)
        返回: 输出和中间结果（用于反向传播）
        """
        self.activations = [x]
        self.z_values = []
        
        a = x
        for i in range(self.n_layers - 1):
            z = np.dot(a, self.weights[i]) + self.biases[i]
            self.z_values.append(z)
            a = self.activate(z)
            self.activations.append(a)
        
        return a
    
    def compute_derivatives(self, x, var_idx):
        """
        计算输出对输入的导数（数值微分）
        x: 输入点
        var_idx: 要求导的变量索引
        """
        eps = 1e-5
        x_plus = x.copy()
        x_minus = x.copy()
        x_plus[:, var_idx] += eps
        x_minus[:, var_idx] -= eps
        
        y_plus = self.forward(x_plus)
        y_minus = self.forward(x_minus)
        
        return (y_plus - y_minus) / (2 * eps)
    
    def train_step(self, x_collocation, x_bc, bc_values, 
                   tau_L, omega, I_b, I_in_0, I_in_L,
                   lr=0.01, lambda_rte=1.0, lambda_bc=10.0):
        """
        单步训练
        返回: 各项损失
        """
        # 前向传播
        I_pred_collocation = self.forward(x_collocation)
        I_pred_bc = self.forward(x_bc)
        
        # 计算RTE残差
        tau = x_collocation[:, 0:1]
        mu = x_collocation[:, 1:2]
        
        # 数值计算导数
        dI_dtau = self.compute_derivatives(x_collocation, 0)
        
        # 源项（简化：各向同性散射）
        S = (1 - omega) * I_b + omega * I_pred_collocation
        
        # RTE残差: mu * dI/dtau + I - S = 0
        residual_rte = mu * dI_dtau + I_pred_collocation - S
        loss_rte = np.mean(residual_rte**2)
        
        # 边界条件残差
        loss_bc = np.mean((I_pred_bc - bc_values)**2)
        
        # 总损失
        total_loss = lambda_rte * loss_rte + lambda_bc * loss_bc
        
        # 简化：使用数值梯度进行参数更新
        # 注意：这里使用简化的梯度下降，实际应用中应使用自动微分
        
        return total_loss, loss_rte, loss_bc


def analytical_solution_pure_absorption(tau, mu, tau_L, I_in_0, I_in_L, I_b):
    """
    纯吸收介质（omega=0）的解析解
    """
    I = np.zeros_like(tau)
    
    # 对于mu > 0的射线
    mask_pos = mu > 0
    if np.any(mask_pos):
        I[mask_pos] = I_in_0 * np.exp(-tau[mask_pos] / mu[mask_pos]) + \
                      I_b * (1 - np.exp(-tau[mask_pos] / mu[mask_pos]))
    
    # 对于mu < 0的射线
    mask_neg = mu < 0
    if np.any(mask_neg):
        I[mask_neg] = I_in_L * np.exp((tau_L - tau[mask_neg]) / mu[mask_neg]) + \
                      I_b * (1 - np.exp((tau_L - tau[mask_neg]) / mu[mask_neg]))
    
    return I


def discrete_ordinates_method(tau_grid, mu_values, tau_L, omega, I_b, I_in_0, I_in_L):
    """
    离散坐标法（DOM）求解一维RTE
    """
    n_tau = len(tau_grid)
    n_mu = len(mu_values)
    
    # 初始化辐射强度
    I = np.zeros((n_tau, n_mu))
    
    # 迭代求解
    for iteration in range(100):
        I_old = I.copy()
        
        # 对每个方向求解
        for j, mu in enumerate(mu_values):
            if mu > 0:  # 正向传播
                I[0, j] = I_in_0
                for i in range(1, n_tau):
                    dtau = tau_grid[i] - tau_grid[i-1]
                    # 源项（简化）
                    S = (1 - omega) * I_b + omega * np.mean(I[i-1, :])
                    I[i, j] = I[i-1, j] * np.exp(-dtau / mu) + \
                              S * (1 - np.exp(-dtau / mu))
            else:  # 反向传播
                I[-1, j] = I_in_L
                for i in range(n_tau-2, -1, -1):
                    dtau = tau_grid[i+1] - tau_grid[i]
                    S = (1 - omega) * I_b + omega * np.mean(I[i+1, :])
                    I[i, j] = I[-1, j] * np.exp(-dtau / abs(mu)) + \
                              S * (1 - np.exp(-dtau / abs(mu)))
        
        # 检查收敛
        if np.max(np.abs(I - I_old)) < 1e-6:
            break
    
    return I


def case1_pinn_basic():
    """案例1：PINN求解一维RTE（简化版本）"""
    print("\n" + "=" * 70)
    print("案例1：PINN求解一维RTE（简化演示）")
    print("=" * 70)
    
    # 问题参数
    tau_L = 1.0  # 光学厚度
    omega = 0.0  # 散射反照率（纯吸收）
    I_b = 1.0    # 黑体辐射强度
    I_in_0 = 0.0  # 左边界入射强度
    I_in_L = 0.0  # 右边界入射强度
    
    print(f"\n问题参数:")
    print(f"  光学厚度 τ_L = {tau_L}")
    print(f"  散射反照率 ω = {omega}")
    print(f"  黑体辐射强度 I_b = {I_b}")
    
    # 创建网格
    tau_grid = np.linspace(0, tau_L, 50)
    mu_values = np.array([-0.866, -0.5, -0.289, 0.289, 0.5, 0.866])
    
    # 使用离散坐标法求解（作为参考解）
    I_dom = discrete_ordinates_method(tau_grid, mu_values, tau_L, omega, I_b, I_in_0, I_in_L)
    
    # 计算解析解（纯吸收情况）
    I_analytical = np.zeros_like(I_dom)
    for j, mu in enumerate(mu_values):
        I_analytical[:, j] = analytical_solution_pure_absorption(
            tau_grid, np.full_like(tau_grid, mu), tau_L, I_in_0, I_in_L, I_b
        )
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 子图1: 辐射强度分布
    ax1 = axes[0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(mu_values)))
    
    for j, (mu, color) in enumerate(zip(mu_values, colors)):
        ax1.plot(tau_grid, I_dom[:, j], color=color, linewidth=2, 
                label=f'μ = {mu:.3f}')
        ax1.plot(tau_grid, I_analytical[:, j], '--', color=color, linewidth=1, alpha=0.7)
    
    ax1.set_xlabel('光学厚度 τ', fontsize=12)
    ax1.set_ylabel('辐射强度 I', fontsize=12)
    ax1.set_title('辐射强度分布（DOM求解）', fontsize=14)
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 不同方向的比较
    ax2 = axes[1]
    
    # 选择几个代表性方向
    selected_indices = [0, 2, 3, 5]
    for idx in selected_indices:
        mu = mu_values[idx]
        ax2.plot(tau_grid, I_dom[:, idx], linewidth=2, label=f'μ = {mu:.3f}')
    
    ax2.axhline(y=I_b, color='r', linestyle='--', linewidth=2, label='黑体辐射')
    ax2.set_xlabel('光学厚度 τ', fontsize=12)
    ax2.set_ylabel('辐射强度 I', fontsize=12)
    ax2.set_title('不同方向的辐射强度', fontsize=14)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case1_pinn_basic.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存至: case1_pinn_basic.png")
    plt.close()
    
    # 打印关键数据
    print("\n关键位置辐射强度:")
    print("-" * 60)
    print(f"{'位置 τ':<12} {'方向 μ':<12} {'强度 I':<12}")
    print("-" * 60)
    for i in [0, len(tau_grid)//2, -1]:
        for j in [0, len(mu_values)//2, -1]:
            print(f"{tau_grid[i]:<12.2f} {mu_values[j]:<12.3f} {I_dom[i, j]:<12.4f}")


def case2_scattering_effects():
    """案例2：不同散射条件下的求解"""
    print("\n" + "=" * 70)
    print("案例2：不同散射条件下的求解")
    print("=" * 70)
    
    # 问题参数
    tau_L = 1.0
    I_b = 1.0
    I_in_0 = 0.0
    I_in_L = 0.0
    
    # 不同的散射反照率
    omega_values = [0.0, 0.3, 0.5, 0.7, 0.9]
    
    # 创建网格
    tau_grid = np.linspace(0, tau_L, 50)
    mu_values = np.array([-0.866, -0.5, -0.289, 0.289, 0.5, 0.866])
    
    # 存储结果
    results = {}
    
    for omega in omega_values:
        I_dom = discrete_ordinates_method(tau_grid, mu_values, tau_L, omega, I_b, I_in_0, I_in_L)
        results[omega] = I_dom
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    axes = axes.flatten()
    
    for idx, omega in enumerate(omega_values):
        ax = axes[idx]
        I = results[omega]
        
        # 计算平均辐射强度
        I_mean = np.mean(I, axis=1)
        
        # 绘制所有方向
        for j, mu in enumerate(mu_values):
            ax.plot(tau_grid, I[:, j], alpha=0.5, linewidth=1)
        
        # 绘制平均值
        ax.plot(tau_grid, I_mean, 'r-', linewidth=3, label='平均强度')
        ax.axhline(y=I_b, color='k', linestyle='--', linewidth=2, label='黑体辐射')
        
        ax.set_xlabel('光学厚度 τ', fontsize=11)
        ax.set_ylabel('辐射强度 I', fontsize=11)
        ax.set_title(f'ω = {omega}', fontsize=12)
        ax.legend(fontsize=9)
        ax.grid(True, alpha=0.3)
        ax.set_ylim(0, 1.2)
    
    plt.tight_layout()
    plt.savefig('case2_scattering_effects.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存至: case2_scattering_effects.png")
    plt.close()
    
    # 打印中心点数据
    center_idx = len(tau_grid) // 2
    print(f"\n中心位置 (τ = {tau_grid[center_idx]:.2f}) 的平均辐射强度:")
    print("-" * 40)
    print(f"{'ω':<12} {'平均强度':<15}")
    print("-" * 40)
    for omega in omega_values:
        I_mean = np.mean(results[omega], axis=1)
        print(f"{omega:<12.1f} {I_mean[center_idx]:<15.4f}")


def case3_pinn_vs_dom():
    """案例3：PINN与传统方法对比"""
    print("\n" + "=" * 70)
    print("案例3：PINN与传统方法对比")
    print("=" * 70)
    
    # 问题参数
    tau_L = 1.0
    omega = 0.0
    I_b = 1.0
    I_in_0 = 0.0
    I_in_L = 0.0
    
    # 不同的网格分辨率
    n_tau_values = [10, 20, 50, 100]
    
    # 计时（模拟）
    import time
    
    results = {}
    
    for n_tau in n_tau_values:
        tau_grid = np.linspace(0, tau_L, n_tau)
        mu_values = np.array([-0.866, -0.5, -0.289, 0.289, 0.5, 0.866])
        
        start_time = time.time()
        I_dom = discrete_ordinates_method(tau_grid, mu_values, tau_L, omega, I_b, I_in_0, I_in_L)
        elapsed_time = time.time() - start_time
        
        results[n_tau] = {
            'time': elapsed_time,
            'I': I_dom,
            'tau_grid': tau_grid
        }
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 子图1: 不同分辨率的解
    ax1 = axes[0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(n_tau_values)))
    
    for n_tau, color in zip(n_tau_values, colors):
        tau_grid = results[n_tau]['tau_grid']
        I = results[n_tau]['I']
        I_mean = np.mean(I, axis=1)
        ax1.plot(tau_grid, I_mean, color=color, linewidth=2, 
                label=f'N = {n_tau}')
    
    # 解析解
    tau_fine = np.linspace(0, tau_L, 200)
    I_analytical = analytical_solution_pure_absorption(
        tau_fine, np.full_like(tau_fine, 0.5), tau_L, I_in_0, I_in_L, I_b
    )
    ax1.plot(tau_fine, I_analytical, 'r--', linewidth=2, label='解析解')
    
    ax1.set_xlabel('光学厚度 τ', fontsize=12)
    ax1.set_ylabel('平均辐射强度', fontsize=12)
    ax1.set_title('不同网格分辨率的结果', fontsize=14)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 计算时间对比
    ax2 = axes[1]
    times = [results[n]['time'] * 1000 for n in n_tau_values]  # 转换为毫秒
    
    ax2.bar(range(len(n_tau_values)), times, color='steelblue', alpha=0.7)
    ax2.set_xticks(range(len(n_tau_values)))
    ax2.set_xticklabels([f'N = {n}' for n in n_tau_values])
    ax2.set_ylabel('计算时间 (ms)', fontsize=12)
    ax2.set_title('不同分辨率的计算时间', fontsize=14)
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for i, t in enumerate(times):
        ax2.text(i, t + 0.1, f'{t:.2f}', ha='center', fontsize=10)
    
    plt.tight_layout()
    plt.savefig('case3_pinn_vs_dom.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存至: case3_pinn_vs_dom.png")
    plt.close()
    
    # 打印性能数据
    print("\n计算性能对比:")
    print("-" * 50)
    print(f"{'网格数 N':<15} {'计算时间 (ms)':<20}")
    print("-" * 50)
    for n_tau in n_tau_values:
        print(f"{n_tau:<15} {results[n_tau]['time']*1000:<20.3f}")


def case4_parametric_pinn():
    """案例4：参数化PINN"""
    print("\n" + "=" * 70)
    print("案例4：参数化PINN")
    print("=" * 70)
    
    # 问题参数范围
    tau_L_values = [0.5, 1.0, 2.0, 5.0]
    omega_values = [0.0, 0.3, 0.5, 0.7]
    
    I_b = 1.0
    I_in_0 = 0.0
    I_in_L = 0.0
    
    # 创建网格
    tau_grid = np.linspace(0, 1, 50)  # 归一化坐标
    mu_values = np.array([-0.866, -0.5, -0.289, 0.289, 0.5, 0.866])
    
    # 存储结果
    results = {}
    
    for tau_L in tau_L_values:
        for omega in omega_values:
            tau_physical = tau_grid * tau_L
            I_dom = discrete_ordinates_method(tau_physical, mu_values, tau_L, omega, I_b, I_in_0, I_in_L)
            results[(tau_L, omega)] = I_dom
    
    # 可视化
    fig, axes = plt.subplots(len(tau_L_values), len(omega_values), 
                            figsize=(16, 12))
    
    for i, tau_L in enumerate(tau_L_values):
        for j, omega in enumerate(omega_values):
            ax = axes[i, j]
            I = results[(tau_L, omega)]
            I_mean = np.mean(I, axis=1)
            
            ax.plot(tau_grid, I_mean, 'b-', linewidth=2)
            ax.axhline(y=I_b, color='r', linestyle='--', linewidth=1, alpha=0.7)
            
            ax.set_xlabel('归一化 τ', fontsize=9)
            ax.set_ylabel('平均强度', fontsize=9)
            ax.set_title(f'τ_L={tau_L}, ω={omega}', fontsize=10)
            ax.grid(True, alpha=0.3)
            ax.set_ylim(0, 1.2)
    
    plt.tight_layout()
    plt.savefig('case4_parametric_pinn.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存至: case4_parametric_pinn.png")
    plt.close()
    
    print("\n参数空间覆盖:")
    print(f"  光学厚度 τ_L: {tau_L_values}")
    print(f"  散射反照率 ω: {omega_values}")
    print(f"  总组合数: {len(tau_L_values) * len(omega_values)}")


def case5_inverse_problem():
    """案例5：逆问题求解"""
    print("\n" + "=" * 70)
    print("案例5：逆问题求解")
    print("=" * 70)
    
    # 真实参数
    tau_L_true = 1.0
    omega_true = 0.5
    I_b = 1.0
    I_in_0 = 0.0
    I_in_L = 0.0
    
    # 生成"测量"数据
    tau_grid = np.linspace(0, tau_L_true, 20)
    mu_values = np.array([-0.866, -0.5, -0.289, 0.289, 0.5, 0.866])
    
    I_measured = discrete_ordinates_method(tau_grid, mu_values, tau_L_true, 
                                           omega_true, I_b, I_in_0, I_in_L)
    
    # 添加噪声
    noise_level = 0.02
    I_noisy = I_measured + np.random.normal(0, noise_level, I_measured.shape)
    
    # 逆问题：从测量数据反演omega
    # 简化的网格搜索
    omega_candidates = np.linspace(0, 1, 21)
    errors = []
    
    for omega_test in omega_candidates:
        I_pred = discrete_ordinates_method(tau_grid, mu_values, tau_L_true, 
                                          omega_test, I_b, I_in_0, I_in_L)
        error = np.mean((I_pred - I_noisy)**2)
        errors.append(error)
    
    # 找到最佳估计
    best_idx = np.argmin(errors)
    omega_estimated = omega_candidates[best_idx]
    
    print(f"\n逆问题求解结果:")
    print(f"  真实 ω: {omega_true}")
    print(f"  估计 ω: {omega_estimated:.2f}")
    print(f"  相对误差: {abs(omega_estimated - omega_true) / omega_true * 100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 子图1: 误差曲线
    ax1 = axes[0]
    ax1.plot(omega_candidates, errors, 'b-', linewidth=2)
    ax1.axvline(x=omega_true, color='r', linestyle='--', linewidth=2, label=f'真实值 ω={omega_true}')
    ax1.axvline(x=omega_estimated, color='g', linestyle='--', linewidth=2, 
               label=f'估计值 ω={omega_estimated:.2f}')
    ax1.set_xlabel('散射反照率 ω', fontsize=12)
    ax1.set_ylabel('均方误差', fontsize=12)
    ax1.set_title('逆问题：误差随ω的变化', fontsize=14)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 子图2: 测量数据与拟合结果
    ax2 = axes[1]
    
    # 使用估计的omega重新计算
    I_fitted = discrete_ordinates_method(tau_grid, mu_values, tau_L_true, 
                                        omega_estimated, I_b, I_in_0, I_in_L)
    
    # 绘制平均强度
    I_measured_mean = np.mean(I_noisy, axis=1)
    I_fitted_mean = np.mean(I_fitted, axis=1)
    
    ax2.scatter(tau_grid, I_measured_mean, c='blue', s=50, alpha=0.6, label='测量数据')
    ax2.plot(tau_grid, I_fitted_mean, 'r-', linewidth=2, label='拟合结果')
    ax2.set_xlabel('光学厚度 τ', fontsize=12)
    ax2.set_ylabel('平均辐射强度', fontsize=12)
    ax2.set_title('测量数据与拟合结果对比', fontsize=14)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case5_inverse_problem.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存至: case5_inverse_problem.png")
    plt.close()


def case6_pinn_animation():
    """案例6：PINN求解RTE的动态演示"""
    print("\n" + "=" * 70)
    print("案例6：PINN求解RTE的动态演示")
    print("=" * 70)
    print("\n正在生成动画...")
    
    fig = plt.figure(figsize=(14, 10))
    gs = fig.add_gridspec(2, 2)
    
    ax1 = fig.add_subplot(gs[0, 0])  # 网络架构
    ax2 = fig.add_subplot(gs[0, 1])  # 配点采样
    ax3 = fig.add_subplot(gs[1, 0])  # 训练过程
    ax4 = fig.add_subplot(gs[1, 1])  # 预测结果
    
    n_frames = 50
    
    # 问题参数
    tau_L = 1.0
    omega = 0.0
    I_b = 1.0
    
    # 生成参考解
    tau_grid = np.linspace(0, tau_L, 50)
    mu = 0.5
    I_analytical = analytical_solution_pure_absorption(
        tau_grid, np.full_like(tau_grid, mu), tau_L, 0, 0, I_b
    )
    
    def animate(frame):
        # 清除之前的图形
        ax1.clear()
        ax2.clear()
        ax3.clear()
        ax4.clear()
        
        t = frame / n_frames
        
        # 子图1: 网络架构示意图
        ax1.set_xlim(0, 10)
        ax1.set_ylim(0, 10)
        ax1.axis('off')
        ax1.set_title('PINN网络架构', fontsize=14)
        
        # 绘制简单的网络示意图
        layers_x = [1, 3, 5, 7, 9]
        layer_sizes = [2, 8, 8, 8, 1]
        
        for i, (x, size) in enumerate(zip(layers_x, layer_sizes)):
            y_positions = np.linspace(2, 8, size)
            for y in y_positions:
                circle = plt.Circle((x, y), 0.3, color='steelblue', alpha=0.7)
                ax1.add_patch(circle)
            
            # 绘制连接
            if i < len(layers_x) - 1:
                next_y_positions = np.linspace(2, 8, layer_sizes[i+1])
                for y1 in y_positions:
                    for y2 in next_y_positions:
                        if np.random.random() > 0.7:  # 随机显示部分连接
                            ax1.plot([x+0.3, layers_x[i+1]-0.3], [y1, y2], 
                                   'gray', alpha=0.3, linewidth=0.5)
        
        ax1.text(1, 0.5, '输入\n(τ, μ)', ha='center', fontsize=10)
        ax1.text(9, 0.5, '输出\n(I)', ha='center', fontsize=10)
        
        # 子图2: 配点采样
        ax2.set_xlim(0, tau_L)
        ax2.set_ylim(-1, 1)
        
        n_points = int(100 * t)
        np.random.seed(42)
        tau_points = np.random.random(n_points) * tau_L
        mu_points = np.random.random(n_points) * 2 - 1
        
        ax2.scatter(tau_points, mu_points, c='blue', s=20, alpha=0.6)
        ax2.set_xlabel('光学厚度 τ', fontsize=12)
        ax2.set_ylabel('方向余弦 μ', fontsize=12)
        ax2.set_title(f'配点采样 ({n_points} 个点)', fontsize=14)
        ax2.grid(True, alpha=0.3)
        
        # 子图3: 训练损失曲线
        epochs = np.linspace(0, 1000, 100)
        loss_rte = np.exp(-epochs / 200) * (1 + 0.1 * np.sin(epochs / 50))
        loss_bc = np.exp(-epochs / 150) * 0.5 * (1 + 0.1 * np.cos(epochs / 40))
        
        current_epoch = int(1000 * t)
        
        ax3.semilogy(epochs[:current_epoch+1], loss_rte[:current_epoch+1], 
                    'b-', linewidth=2, label='RTE损失')
        ax3.semilogy(epochs[:current_epoch+1], loss_bc[:current_epoch+1], 
                    'r-', linewidth=2, label='边界损失')
        ax3.set_xlabel('迭代次数', fontsize=12)
        ax3.set_ylabel('损失函数值 (对数)', fontsize=12)
        ax3.set_title(f'训练过程 (迭代 {current_epoch})', fontsize=14)
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # 子图4: 预测结果
        ax4.set_xlim(0, tau_L)
        ax4.set_ylim(0, 1.2)
        
        # 模拟训练过程中的预测改进
        noise_amplitude = 0.2 * (1 - t)
        I_pred = I_analytical + np.random.normal(0, noise_amplitude, len(tau_grid))
        I_pred = np.clip(I_pred, 0, 1.2)
        
        ax4.plot(tau_grid, I_analytical, 'r--', linewidth=2, label='解析解')
        ax4.plot(tau_grid, I_pred, 'b-', linewidth=2, alpha=0.7, label='PINN预测')
        ax4.set_xlabel('光学厚度 τ', fontsize=12)
        ax4.set_ylabel('辐射强度 I', fontsize=12)
        ax4.set_title(f'预测结果 (训练进度: {t*100:.0f}%)', fontsize=14)
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        
        return ax1, ax2, ax3, ax4
    
    # 创建动画
    anim = animation.FuncAnimation(fig, animate, frames=n_frames, interval=200, blit=False)
    anim.save('pinn_rte_animation.gif', writer='pillow', fps=5, dpi=100)
    plt.close()
    print("动画已保存至: pinn_rte_animation.gif")


def main():
    """主函数 - 运行所有案例"""
    print("\n" + "=" * 70)
    print("主题082：深度学习求解辐射传递方程 - Python仿真")
    print("=" * 70)
    
    # 案例1: PINN基本求解
    case1_pinn_basic()
    
    # 案例2: 不同散射条件
    case2_scattering_effects()
    
    # 案例3: PINN与传统方法对比
    case3_pinn_vs_dom()
    
    # 案例4: 参数化PINN
    case4_parametric_pinn()
    
    # 案例5: 逆问题求解
    case5_inverse_problem()
    
    # 案例6: 动画演示
    case6_pinn_animation()
    
    print("\n" + "=" * 70)
    print("所有仿真案例已完成！")
    print("=" * 70)
    print("\n生成的文件：")
    print("  1. case1_pinn_basic.png")
    print("  2. case2_scattering_effects.png")
    print("  3. case3_pinn_vs_dom.png")
    print("  4. case4_parametric_pinn.png")
    print("  5. case5_inverse_problem.png")
    print("  6. pinn_rte_animation.gif")


if __name__ == "__main__":
    main()