下载地址：http://lanzou.com.cn/i60ca2287

📁 output/jisuanmoxingaiban/
├── 📄 README.md188 B
├── 📄 pom.xml1.4 KB
├── 📄 package.json688 B
├── 📄 config/application.properties603 B
├── 📄 infer/Factory.php3.7 KB
├── 📄 infer/Provider.sql3.1 KB
├── 📄 weights/Engine.js3.9 KB
├── 📄 weights/Validator.ts3.6 KB
├── 📄 view/Controller.go3.1 KB
├── 📄 infer/Registry.java4.2 KB
├── 📄 src/main/java/Proxy.java7 KB
├── 📄 utils/Scheduler.cpp1.6 KB
├── 📄 src/main/java/Executor.java3.9 KB
├── 📄 utils/Buffer.py5.6 KB
├── 📄 view/Server.py5.9 KB
├── 📄 view/Manager.go4.1 KB
├── 📄 weights/Parser.js4.5 KB
├── 📄 utils/Client.js4.7 KB
├── 📄 view/Transformer.cpp1.5 KB
├── 📄 config/Worker.json688 B
├── 📄 utils/Util.py5.6 KB
├── 📄 src/main/java/Helper.java6.7 KB
├── 📄 view/Repository.php2.8 KB
├── 📄 config/Dispatcher.properties604 B
├── 📄 config/Loader.xml1.4 KB

项目编译入口：

Project Structure

Project : 仿真计算模型AI版

Folder : jisuanmoxingaiban

Files : 26

Size : 80.9 KB

Generated: 2026-03-23 18:47:36

jisuanmoxingaiban/
├── README.md [188 B]
├── config/
│ ├── Dispatcher.properties [604 B]
│ ├── Loader.xml [1.4 KB]
│ ├── Worker.json [688 B]
│ └── application.properties [603 B]
├── factory/
├── infer/
│ ├── Factory.php [3.7 KB]
│ ├── Provider.sql [3.1 KB]
│ └── Registry.java [4.2 KB]
├── lib/
├── package.json [688 B]
├── pom.xml [1.4 KB]
├── src/
│ ├── main/
│ │ ├── java/
│ │ │ ├── Executor.java [3.9 KB]
│ │ │ ├── Helper.java [6.7 KB]
│ │ │ └── Proxy.java [7 KB]
│ │ └── resources/
│ └── test/
│ └── java/
├── utils/
│ ├── Buffer.py [5.6 KB]
│ ├── Client.js [4.7 KB]
│ ├── Scheduler.cpp [1.6 KB]
│ └── Util.py [5.6 KB]
├── view/
│ ├── Controller.go [3.1 KB]
│ ├── Manager.go [4.1 KB]
│ ├── Repository.php [2.8 KB]
│ ├── Server.py [5.9 KB]
│ └── Transformer.cpp [1.5 KB]
└── weights/
├── Engine.js [3.9 KB]
├── Parser.js [4.5 KB]
└── Validator.ts [3.6 KB]

python
"""
thermal_simulation_base.py
基础热传导仿真模块，包含有限差分求解器、边界条件处理与可视化
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from typing import Tuple, Callable, Optional, Dict, Any
import dataclasses
from scipy.sparse import diags, csr_matrix
from scipy.sparse.linalg import spsolve
import logging
import time

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(name)

@dataclasses.dataclass
class SimulationConfig:
"""仿真配置数据类"""
nx: int = 100 # 空间网格数
ny: int = 100 # y方向网格数（2D情况）
nt: int = 500 # 时间步数
Lx: float = 1.0 # x方向长度
Ly: float = 1.0 # y方向长度
alpha: float = 0.01 # 热扩散系数
dt: float = None # 时间步长（自动计算稳定性条件）
method: str = "explicit" # explicit, implicit, crank_nicolson
dim: int = 2 # 维度: 1 或 2

class ThermalSolver2D:
"""
二维热传导方程求解器
支持显式、隐式和Crank-Nicolson格式
"""

def __init__(self, config: SimulationConfig):
    self.config = config
    self.dx = config.Lx / (config.nx - 1)
    self.dy = config.Ly / (config.ny - 1)

    # 稳定性条件自动确定时间步长
    if config.dt is None:
        dx2 = self.dx ** 2
        dy2 = self.dy ** 2
        # 显式格式稳定性条件
        max_dt = 0.5 * min(dx2, dy2) / config.alpha
        self.dt = 0.8 * max_dt  # 安全系数
        logger.info(f"自动设置时间步长 dt = {self.dt:.6f}")
    else:
        self.dt = config.dt

    # 创建网格坐标
    self.x = np.linspace(0, config.Lx, config.nx)
    self.y = np.linspace(0, config.Ly, config.ny)
    self.X, self.Y = np.meshgrid(self.x, self.y)

    # 初始化温度场
    self.u = np.zeros((config.nx, config.ny))
    self.u_prev = np.zeros((config.nx, config.ny))

    # 缓存隐式格式的矩阵（避免重复构建）
    self._A_matrix = None
    self._B_matrix = None

def set_initial_condition(self, func: Callable[[np.ndarray, np.ndarray], np.ndarray]) -> None:
    """设置初始温度分布"""
    self.u = func(self.X, self.Y)
    self.u_prev = self.u.copy()

def set_boundary_conditions(self, bc_type: str = "dirichlet", 
                             bc_values: Dict[str, float] = None) -> None:
    """
    设置边界条件
    bc_type: "dirichlet", "neumann", "mixed"
    bc_values: 边界值字典，如 {'left': 0, 'right': 0, 'top': 0, 'bottom': 0}
    """
    self.bc_type = bc_type
    self.bc_values = bc_values or {}

def _build_implicit_matrices(self):
    """构建隐式和Crank-Nicolson格式所需的稀疏矩阵"""
    nx, ny = self.config.nx, self.config.ny
    N = nx * ny

    dx2 = self.dx ** 2
    dy2 = self.dy ** 2

    # 有限差分系数
    rx = self.config.alpha * self.dt / dx2
    ry = self.config.alpha * self.dt / dy2

    # 构建五点差分格式的稀疏矩阵
    diagonals = []
    offsets = []

    # 主对角线
    main_diag = 1 + 2 * rx + 2 * ry
    diagonals.append(main_diag * np.ones(N))
    offsets.append(0)

    # x方向相邻点
    diagonals.append(-rx * np.ones(N - 1))
    offsets.append(1)
    diagonals.append(-rx * np.ones(N - 1))
    offsets.append(-1)

    # y方向相邻点
    diagonals.append(-ry * np.ones(N - nx))
    offsets.append(nx)
    diagonals.append(-ry * np.ones(N - nx))
    offsets.append(-nx)

    self._A_matrix = diags(diagonals, offsets, format='csr')

    # 对于Crank-Nicolson，需要右端矩阵
    if self.config.method == "crank_nicolson":
        diagonals_b = []
        main_diag_b = 1 - 2 * rx - 2 * ry
        diagonals_b.append(main_diag_b * np.ones(N))
        diagonals_b.append(rx * np.ones(N - 1))
        diagonals_b.append(rx * np.ones(N - 1))
        diagonals_b.append(ry * np.ones(N - nx))
        diagonals_b.append(ry * np.ones(N - nx))
        self._B_matrix = diags(diagonals_b, offsets, format='csr')

def _apply_boundary_conditions_to_matrix(self, matrix: csr_matrix) -> csr_matrix:
    """将边界条件应用到矩阵"""
    nx, ny = self.config.nx, self.config.ny
    matrix = matrix.tolil()

    # 获取边界索引
    left_boundary = [i * nx for i in range(ny)]
    right_boundary = [i * nx + (nx - 1) for i in range(ny)]
    bottom_boundary = list(range(nx))
    top_boundary = [nx * (ny - 1) + i for i in range(nx)]

    all_boundary = left_boundary + right_boundary + bottom_boundary + top_boundary

    for idx in all_boundary:
        if self.bc_type == "dirichlet":
            # Dirichlet: 设置对角线为1，其他为0
            matrix[idx, :] = 0
            matrix[idx, idx] = 1
        elif self.bc_type == "neumann":
            # Neumann: 一阶近似，待简化
            pass

    return matrix.tocsr()

def step_explicit(self) -> np.ndarray:
    """显式时间推进"""
    nx, ny = self.config.nx, self.config.ny
    dx2 = self.dx ** 2
    dy2 = self.dy ** 2

    u_new = self.u.copy()

    # 内部点更新
    u_new[1:-1, 1:-1] = self.u[1:-1, 1:-1] + \
        self.config.alpha * self.dt * (
            (self.u[2:, 1:-1] - 2 * self.u[1:-1, 1:-1] + self.u[:-2, 1:-1]) / dx2 +
            (self.u[1:-1, 2:] - 2 * self.u[1:-1, 1:-1] + self.u[1:-1, :-2]) / dy2
        )

    # 边界条件处理
    self._apply_boundary_conditions(u_new)

    self.u = u_new
    return self.u

def _apply_boundary_conditions(self, u_field: np.ndarray) -> None:
    """应用边界条件到温度场"""
    if self.bc_type == "dirichlet":
        if 'left' in self.bc_values:
            u_field[0, :] = self.bc_values['left']
        if 'right' in self.bc_values:
            u_field[-1, :] = self.bc_values['right']
        if 'bottom' in self.bc_values:
            u_field[:, 0] = self.bc_values['bottom']
        if 'top' in self.bc_values:
            u_field[:, -1] = self.bc_values['top']
    elif self.bc_type == "neumann":
        # 简化处理
        if 'left' in self.bc_values:
            u_field[0, :] = u_field[1, :] - self.bc_values['left'] * self.dx

def step_implicit(self) -> np.ndarray:
    """隐式时间推进"""
    if self._A_matrix is None:
        self._build_implicit_matrices()

    nx, ny = self.config.nx, self.config.ny
    N = nx * ny

    # 将温度场展平为向量
    u_flat = self.u.flatten()

    # 构建右端项（隐式格式中就是当前时刻的u）
    rhs = u_flat.copy()

    # 应用边界条件
    A = self._apply_boundary_conditions_to_matrix(self._A_matrix)

    # 求解线性方程组
    u_new_flat = spsolve(A, rhs)
    self.u = u_new_flat.reshape((nx, ny))

    return self.u

def step_crank_nicolson(self) -> np.ndarray:
    """Crank-Nicolson时间推进"""
    if self._A_matrix is None or self._B_matrix is None:
        self._build_implicit_matrices()

    nx, ny = self.config.nx, self.config.ny
    N = nx * ny

    u_flat = self.u.flatten()

    # 右端项 = B * u_current
    rhs = self._B_matrix.dot(u_flat)

    # 求解 A * u_new = rhs
    A = self._apply_boundary_conditions_to_matrix(self._A_matrix)
    u_new_flat = spsolve(A, rhs)
    self.u = u_new_flat.reshape((nx, ny))

    return self.u

def step(self) -> np.ndarray:
    """执行一个时间步"""
    if self.config.method == "explicit":
        return self.step_explicit()
    elif self.config.method == "implicit":
        return self.step_implicit()
    elif self.config.method == "crank_nicolson":
        return self.step_crank_nicolson()
    else:
        raise ValueError(f"未知方法: {self.config.method}")

def run(self, callback: Optional[Callable] = None) -> List[np.ndarray]:
    """运行完整仿真，可选回调函数记录中间状态"""
    history = []
    start_time = time.time()

    for step_idx in range(self.config.nt):
        self.step()

        if callback is not None:
            callback(step_idx, self.u)

        # 每100步记录一次历史（节省内存）
        if step_idx % 100 == 0:
            history.append(self.u.copy())
            logger.info(f"步数 {step_idx}/{self.config.nt}, 最大温度: {np.max(self.u):.4f}")

    elapsed = time.time() - start_time
    logger.info(f"仿真完成，耗时 {elapsed:.2f} 秒")

    return history

def animate(self, history: List[np.ndarray], interval: int = 50) -> FuncAnimation:
    """动画显示温度场演化"""
    fig, ax = plt.subplots(figsize=(8, 6))

    def update(frame):
        ax.clear()
        im = ax.contourf(self.X, self.Y, history[frame].T, levels=50, cmap='hot')
        ax.set_title(f'温度分布 - 时间步 {frame * 100}')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        return im,

    anim = FuncAnimation(fig, update, frames=len(history), interval=interval)
    plt.colorbar(ax.contourf(self.X, self.Y, history[0].T, levels=50, cmap='hot'), ax=ax)
    return anim

示例使用

if name == "main":

# 创建配置
config = SimulationConfig(
    nx=80, ny=80, nt=1000, method="explicit"
)

solver = ThermalSolver2D(config)

# 初始条件：中心高温区域
def initial_temp(x, y):
    return np.exp(-((x-0.5)**2 + (y-0.5)**2) / 0.05)

solver.set_initial_condition(initial_temp)
solver.set_boundary_conditions("dirichlet", {'left': 0, 'right': 0, 'bottom': 0, 'top': 0})

# 运行仿真
history = solver.run()

# 显示最终结果
plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.contourf(solver.X, solver.Y, history[0].T, levels=50, cmap='hot')
plt.title('初始温度')
plt.colorbar()

plt.subplot(122)
plt.contourf(solver.X, solver.Y, history[-1].T, levels=50, cmap='hot')
plt.title(f'最终温度 (t={config.nt * solver.dt:.2f})')
plt.colorbar()
plt.tight_layout()
plt.show()

第二章 AI增强：物理信息神经网络代理模型
传统数值仿真在精细网格下计算量巨大。本节引入物理信息神经网络作为代理模型，在保留物理约束的同时实现快速预测。

python
"""
pinn_thermal_surrogate.py
基于物理信息神经网络的热传导代理模型
"""

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
import numpy as np
import matplotlib.pyplot as plt
from typing import Tuple, List, Callable, Optional
import itertools
from tqdm import tqdm

class PhysicsInformedNN(nn.Module):
"""
物理信息神经网络
同时拟合数据并满足热传导方程残差
"""

def __init__(self, layers: List[int], alpha: float = 0.01):
    """
    参数:
        layers: 网络结构，如 [2, 50, 50, 50, 1]
        alpha: 热扩散系数
    """
    super().__init__()
    self.alpha = alpha

    # 构建全连接网络
    self.network = nn.Sequential()
    for i in range(len(layers) - 2):
        self.network.append(nn.Linear(layers[i], layers[i+1]))
        self.network.append(nn.Tanh())
    self.network.append(nn.Linear(layers[-2], layers[-1]))

    # 权重初始化
    self._initialize_weights()

def _initialize_weights(self):
    for m in self.modules():
        if isinstance(m, nn.Linear):
            nn.init.xavier_normal_(m.weight)
            nn.init.zeros_(m.bias)

def forward(self, x: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
    """前向传播，输入为空间坐标和时间"""
    # 拼接输入 [x, y, t] 或 [x, t]
    if x.dim() == 1:
        x = x.unsqueeze(1)
    if t.dim() == 1:
        t = t.unsqueeze(1)
    inputs = torch.cat([x, t], dim=1)
    return self.network(inputs)

def compute_pde_residual(self, x: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
    """
    计算热传导方程残差: u_t - alpha * u_xx
    使用自动微分
    """
    # 启用梯度
    x.requires_grad_(True)
    t.requires_grad_(True)

    u = self.forward(x, t)

    # 一阶时间导数
    u_t = torch.autograd.grad(
        u, t, grad_outputs=torch.ones_like(u),
        create_graph=True, retain_graph=True
    )[0]

    # 二阶空间导数
    u_x = torch.autograd.grad(
        u, x, grad_outputs=torch.ones_like(u),
        create_graph=True, retain_graph=True
    )[0]
    u_xx = torch.autograd.grad(
        u_x, x, grad_outputs=torch.ones_like(u_x),
        create_graph=True, retain_graph=True
    )[0]

    # PDE残差
    residual = u_t - self.alpha * u_xx
    return residual

class PINNThermalSolver:
"""
使用PINN求解热传导问题的封装类
支持边界条件硬约束和软约束
"""

def __init__(self, layers: List[int], alpha: float, domain: Tuple[float, float, float, float]):
    """
    参数:
        layers: 网络结构
        alpha: 热扩散系数
        domain: (x_min, x_max, t_min, t_max)
    """
    self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    self.model = PhysicsInformedNN(layers, alpha).to(self.device)
    self.alpha = alpha
    self.x_min, self.x_max, self.t_min, self.t_max = domain

def generate_collocation_points(self, n_x: int, n_t: int) -> Tuple[torch.Tensor, torch.Tensor]:
    """生成配置点用于PDE残差计算"""
    x = torch.linspace(self.x_min, self.x_max, n_x, device=self.device)
    t = torch.linspace(self.t_min, self.t_max, n_t, device=self.device)

    # 网格化
    X, T = torch.meshgrid(x, t, indexing='ij')
    X_flat = X.reshape(-1, 1)
    T_flat = T.reshape(-1, 1)

    return X_flat, T_flat

def generate_boundary_points(self, n_points: int) -> Dict[str, Tuple[torch.Tensor, torch.Tensor]]:
    """生成边界点"""
    points = {}

    # 左边界 x = x_min
    t_left = torch.linspace(self.t_min, self.t_max, n_points, device=self.device)
    x_left = torch.full_like(t_left, self.x_min)
    points['left'] = (x_left.unsqueeze(1), t_left.unsqueeze(1))

    # 右边界 x = x_max
    t_right = torch.linspace(self.t_min, self.t_max, n_points, device=self.device)
    x_right = torch.full_like(t_right, self.x_max)
    points['right'] = (x_right.unsqueeze(1), t_right.unsqueeze(1))

    # 初始条件 t = t_min
    x_initial = torch.linspace(self.x_min, self.x_max, n_points, device=self.device)
    t_initial = torch.full_like(x_initial, self.t_min)
    points['initial'] = (x_initial.unsqueeze(1), t_initial.unsqueeze(1))

    return points

def train(self, 
          n_collocation: int = 10000,
          n_boundary: int = 100,
          epochs: int = 10000,
          lr: float = 1e-3,
          pde_weight: float = 1.0,
          bc_weight: float = 1.0,
          ic_weight: float = 1.0,
          verbose: bool = True) -> List[float]:
    """
    训练PINN模型
    损失函数 = pde_weight * L_pde + bc_weight * L_bc + ic_weight * L_ic
    """
    # 生成配置点和边界点
    n_x_coll = int(np.sqrt(n_collocation))
    n_t_coll = n_x_coll
    X_coll, T_coll = self.generate_collocation_points(n_x_coll, n_t_coll)

    boundary_points = self.generate_boundary_points(n_boundary)

    # 初始条件的目标值（需要根据具体问题定义）
    # 这里使用高斯脉冲作为示例
    def initial_condition_func(x):
        return torch.exp(-((x - 0.5)**2) / 0.05)

    # 边界条件目标值
    def boundary_condition_func(x, t):
        return torch.zeros_like(x)

    # 优化器
    optimizer = optim.Adam(self.model.parameters(), lr=lr)
    scheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience=500, factor=0.5)

    losses = []

    progress_bar = tqdm(range(epochs), desc="Training PINN", disable=not verbose)

    for epoch in progress_bar:
        optimizer.zero_grad()

        # 1. PDE残差损失
        X_coll.requires_grad_(True)
        T_coll.requires_grad_(True)
        residual = self.model.compute_pde_residual(X_coll, T_coll)
        loss_pde = torch.mean(residual ** 2)

        # 2. 边界条件损失
        loss_bc = 0.0
        for boundary_name, (x_b, t_b) in boundary_points.items():
            if boundary_name != 'initial':  # 初始条件单独处理
                u_pred = self.model(x_b, t_b)
                u_target = boundary_condition_func(x_b, t_b)
                loss_bc += torch.mean((u_pred - u_target) ** 2)

        # 3. 初始条件损失
        x_ic, t_ic = boundary_points['initial']
        u_ic_pred = self.model(x_ic, t_ic)
        u_ic_target = initial_condition_func(x_ic)
        loss_ic = torch.mean((u_ic_pred - u_ic_target) ** 2)

        # 总损失
        total_loss = pde_weight * loss_pde + bc_weight * loss_bc + ic_weight * loss_ic
        total_loss.backward()

        optimizer.step()
        scheduler.step(total_loss)

        losses.append(total_loss.item())

        if verbose and epoch % 500 == 0:
            progress_bar.set_postfix({
                'loss': f'{total_loss.item():.4e}',
                'pde': f'{loss_pde.item():.4e}',
                'bc': f'{loss_bc.item():.4e}'
            })

    return losses

def predict(self, x: np.ndarray, t: np.ndarray) -> np.ndarray:
    """在给定点进行预测"""
    self.model.eval()
    with torch.no_grad():
        x_tensor = torch.tensor(x, dtype=torch.float32, device=self.device)
        t_tensor = torch.tensor(t, dtype=torch.float32, device=self.device)
        if x_tensor.dim() == 1:
            x_tensor = x_tensor.unsqueeze(1)
        if t_tensor.dim() == 1:
            t_tensor = t_tensor.unsqueeze(1)
        u_pred = self.model(x_tensor, t_tensor)
        return u_pred.cpu().numpy().flatten()

def solve_pde(self, x_grid: np.ndarray, t_grid: np.ndarray) -> np.ndarray:
    """在网格上求解PDE"""
    X, T = np.meshgrid(x_grid, t_grid, indexing='ij')
    X_flat = X.flatten()
    T_flat = T.flatten()
    U_flat = self.predict(X_flat, T_flat)
    return U_flat.reshape(X.shape)

def visualize_solution(self, x_grid: np.ndarray, t_grid: np.ndarray, 
                       title: str = "PINN求解热传导方程"):
    """可视化求解结果"""
    U = self.solve_pde(x_grid, t_grid)

    fig, axes = plt.subplots(1, 3, figsize=(15, 4))

    # 温度分布云图
    X, T = np.meshgrid(x_grid, t_grid, indexing='ij')
    im = axes[0].contourf(X, T, U, levels=50, cmap='hot')
    axes[0].set_xlabel('x')
    axes[0].set_ylabel('t')
    axes[0].set_title('温度分布 u(x,t)')
    plt.colorbar(im, ax=axes[0])

    # 不同时间截面的温度曲线
    time_indices = [0, len(t_grid)//4, len(t_grid)//2, len(t_grid)-1]
    for idx in time_indices:
        axes[1].plot(x_grid, U[:, idx], label=f't={t_grid[idx]:.2f}')
    axes[1].set_xlabel('x')
    axes[1].set_ylabel('温度')
    axes[1].legend()
    axes[1].set_title('不同时刻温度分布')

    # 不同位置的时间演化
    pos_indices = [0, len(x_grid)//4, len(x_grid)//2, len(x_grid)-1]
    for idx in pos_indices:
        axes[2].plot(t_grid, U[idx, :], label=f'x={x_grid[idx]:.2f}')
    axes[2].set_xlabel('t')
    axes[2].set_ylabel('温度')
    axes[2].legend()
    axes[2].set_title('不同位置温度演化')

    plt.suptitle(title)
    plt.tight_layout()
    plt.show()

高级PINN变体：带自适应权重的多任务学习

class AdaptiveWeightedPINN(PhysicsInformedNN):
"""
使用自适应权重机制的多任务PINN
基于不确定性加权或动态加权
"""

def __init__(self, layers: List[int], alpha: float, 
             adaptive_method: str = "uncertainty"):
    super().__init__(layers, alpha)
    self.adaptive_method = adaptive_method

    # 可学习的损失权重（对数形式保证正性）
    self.log_pde_weight = nn.Parameter(torch.tensor(0.0))
    self.log_bc_weight = nn.Parameter(torch.tensor(0.0))
    self.log_ic_weight = nn.Parameter(torch.tensor(0.0))

def get_weights(self):
    """获取当前权重"""
    return (torch.exp(self.log_pde_weight),
            torch.exp(self.log_bc_weight),
            torch.exp(self.log_ic_weight))

def compute_weighted_loss(self, loss_pde, loss_bc, loss_ic):
    """计算加权损失（自适应）"""
    if self.adaptive_method == "uncertainty":
        # 不确定性加权：L = exp(-log_sigma) * loss + log_sigma
        w_pde = torch.exp(-self.log_pde_weight)
        w_bc = torch.exp(-self.log_bc_weight)
        w_ic = torch.exp(-self.log_ic_weight)

        total_loss = (w_pde * loss_pde + self.log_pde_weight +
                     w_bc * loss_bc + self.log_bc_weight +
                     w_ic * loss_ic + self.log_ic_weight)
    else:
        # 固定权重
        w_pde, w_bc, w_ic = 1.0, 1.0, 1.0
        total_loss = w_pde * loss_pde + w_bc * loss_bc + w_ic * loss_ic

    return total_loss

案例：与数值解对比

def compare_with_numerical():
"""对比PINN解与有限差分数值解"""

# 参数设置
alpha = 0.01
x_min, x_max = 0.0, 1.0
t_max = 0.5

# 1. 使用PINN求解
layers = [2, 50, 50, 50, 1]
pinn_solver = PINNThermalSolver(layers, alpha, (x_min, x_max, 0.0, t_max))

# 训练
pinn_solver.train(n_collocation=5000, n_boundary=100, epochs=3000, verbose=False)

# 2. 使用有限差分求解
from thermal_simulation_base import SimulationConfig, ThermalSolver2D

# 简化为一维问题，使用2D求解器但y方向单网格
config = SimulationConfig(
    nx=100, ny=2, nt=500, 
    Lx=1.0, Ly=0.01,
    alpha=alpha,
    method="explicit"
)

solver = ThermalSolver2D(config)

def initial_1d(x, y):
    return np.exp(-((x-0.5)**2) / 0.05)

solver.set_initial_condition(initial_1d)
solver.set_boundary_conditions("dirichlet", {'left': 0, 'right': 0})

# 运行数值仿真
history = solver.run()

# 提取指定时间的数值解
t_numerical = np.linspace(0, t_max, len(history))
x_grid = np.linspace(0, 1, 100)

# PINN预测
t_pinn = np.linspace(0, t_max, 10)
U_pinn = pinn_solver.solve_pde(x_grid, t_pinn)

# 可视化对比
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

for idx, t_idx in enumerate([0, 2, 5, 8]):
    ax = axes[idx // 2, idx % 2]
    t_val = t_pinn[t_idx]

    # 数值解插值到相同时间
    t_step = int(t_val / solver.dt)
    if t_step < len(history):
        u_numerical = history[t_step][:, 0]  # y方向取第一列
    else:
        u_numerical = history[-1][:, 0]

    ax.plot(x_grid, U_pinn[:, t_idx], 'b-', label='PINN', linewidth=2)
    ax.plot(x_grid, u_numerical, 'r--', label='FDM', linewidth=2)
    ax.set_title(f't = {t_val:.3f}')
    ax.set_xlabel('x')
    ax.set_ylabel('温度')
    ax.legend()
    ax.grid(True, alpha=0.3)

plt.suptitle('PINN与有限差分法对比', fontsize=14)
plt.tight_layout()
plt.show()

# 计算相对误差
error = np.abs(U_pinn - u_numerical.reshape(-1, 1)) / (np.abs(u_numerical.reshape(-1, 1)) + 1e-8)
print(f"平均相对误差: {np.mean(error):.4f}")
print(f"最大相对误差: {np.max(error):.4f}")

if name == "main":
compare_with_numerical()
第三章 不确定性量化：高斯过程回归
在仿真中引入不确定性量化，可以评估模型预测的置信区间，对于高风险决策至关重要。

python
"""
gaussian_process_emulator.py
基于高斯过程的仿真代理模型与不确定性量化
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cholesky, solve_triangular
from scipy.spatial.distance import cdist
from sklearn.preprocessing import StandardScaler
from typing import Tuple, List, Optional, Callable, Dict
import warnings
import itertools
from dataclasses import dataclass

@dataclass
class GPConfig:
"""高斯过程配置"""
length_scale: float = 1.0
variance: float = 1.0
noise_variance: float = 1e-6
kernel: str = "rbf" # rbf, matern32, matern52, rational_quadratic

class GaussianProcessEmulator:
"""
高斯过程回归仿真代理模型
提供预测均值、方差和置信区间
"""

def __init__(self, config: GPConfig):
    self.config = config
    self.X_train = None
    self.y_train = None
    self.alpha = None  # 预计算的系数
    self.L = None      # Cholesky分解下三角矩阵
    self.scaler_X = StandardScaler()
    self.scaler_y = StandardScaler()

def _kernel_rbf(self, X1: np.ndarray, X2: np.ndarray) -> np.ndarray:
    """径向基函数核"""
    sqdist = cdist(X1 / self.config.length_scale, 
                   X2 / self.config.length_scale, 
                   metric='sqeuclidean')
    return self.config.variance * np.exp(-0.5 * sqdist)

def _kernel_matern32(self, X1: np.ndarray, X2: np.ndarray) -> np.ndarray:
    """Matérn 3/2 核"""
    dist = cdist(X1 / self.config.length_scale, 
                 X2 / self.config.length_scale, 
                 metric='euclidean')
    sqrt3_dist = np.sqrt(3) * dist
    return self.config.variance * (1 + sqrt3_dist) * np.exp(-sqrt3_dist)

def _kernel_matern52(self, X1: np.ndarray, X2: np.ndarray) -> np.ndarray:
    """Matérn 5/2 核"""
    dist = cdist(X1 / self.config.length_scale, 
                 X2 / self.config.length_scale, 
                 metric='euclidean')
    sqrt5_dist = np.sqrt(5) * dist
    return self.config.variance * (1 + sqrt5_dist + (5/3) * dist**2) * np.exp(-sqrt5_dist)

def _kernel_rational_quadratic(self, X1: np.ndarray, X2: np.ndarray, alpha: float = 1.0) -> np.ndarray:
    """有理二次核"""
    sqdist = cdist(X1 / self.config.length_scale, 
                   X2 / self.config.length_scale, 
                   metric='sqeuclidean')
    return self.config.variance * (1 + sqdist / (2 * alpha)) ** (-alpha)

def compute_kernel(self, X1: np.ndarray, X2: np.ndarray) -> np.ndarray:
    """根据配置选择核函数"""
    if self.config.kernel == "rbf":
        return self._kernel_rbf(X1, X2)
    elif self.config.kernel == "matern32":
        return self._kernel_matern32(X1, X2)
    elif self.config.kernel == "matern52":
        return self._kernel_matern52(X1, X2)
    elif self.config.kernel == "rational_quadratic":
        return self._kernel_rational_quadratic(X1, X2)
    else:
        raise ValueError(f"未知核函数: {self.config.kernel}")

def fit(self, X: np.ndarray, y: np.ndarray) -> 'GaussianProcessEmulator':
    """
    训练高斯过程模型
    X: (n_samples, n_features)
    y: (n_samples,)
    """
    self.X_train = X.copy()

    # 标准化输入和输出
    X_sc