2D Canvas热方程模拟出现异常图案及数值爆炸问题求助

问题根源分析 #

• 原地更新网格导致数值污染：当前代码计算新温度时直接覆盖原网格值，后续单元格计算会使用已经更新的邻域数据，误差快速累积后引发数值爆炸。
• 热源数组Q的边界缺失：Q数组仅初始化到width-2和height-2，边缘单元格（i=width-1或j=height-1）对应的Q值为undefined，参与计算后产生NaN或异常数值。
• 温度值无约束导致颜色渲染异常：当温度超过100时，HSL色调计算会出现负数，引发不符合预期的颜色显示。

修复方案 #

1. 使用临时网格存储新值：计算温度更新时，先将所有新值存入临时数组，计算完成后再替换原网格，避免邻域数据被提前覆盖。
2. 补全Q数组的边界初始化：将Q数组的循环范围改为i < width和j < height，确保所有单元格都有对应的热源值。
3. 约束温度范围（可选但推荐）：对计算后的温度值进行范围限制，避免数值无限增长，同时保证颜色渲染的稳定性。
4. 验证稳定性条件：2D显式热方程的稳定性要求为α*dt*(1/dx² + 1/dy²) ≤ 0.5，当前参数符合要求，但后续调整参数时需确保该条件始终成立。

修复后的完整代码 #

const width = 200; // 网格宽度
const height = 200; // 网格高度
const dt = 0.25; // 时间步长
const dx = 1; // x方向空间步长
const dy = 1; // y方向空间步长
const alpha = 0.25; // 热扩散系数

const Q = [];
// 热源强度
const Q0 = 80;

const r1 = 8;
const r2 = 12;

// 补全Q数组的所有边界单元格
for (let i = 0; i < width; i++) {
    Q[i] = [];
    for (let j = 0; j < height; j++) {
        // 计算到中心的距离
        const r = Math.sqrt((i - width / 2) ** 2 + (j - height / 2) ** 2);
        Q[i][j] = (r1 < r && r < r2) ? Q0 : 0;
    }
}

let grid = []; // 存储温度网格

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');

// 初始化网格温度为50
for (let i = 0; i < width; i++) {
    grid[i] = [];
    for (let j = 0; j < height; j++) {
        grid[i][j] = 50;
    }
}

function updateGrid() {
    // 创建临时网格存储新温度值
    const newGrid = [];
    // 复制原网格的边界值（边界采用绝热条件，温度不变）
    for (let i = 0; i < width; i++) {
        newGrid[i] = [];
        for (let j = 0; j < height; j++) {
            newGrid[i][j] = grid[i][j];
        }
    }

    // 更新内部单元格的温度
    for (let i = 1; i < width - 1; i++) {
        for (let j = 1; j < height - 1; j++) {
            const d2Tdx2 = (grid[i + 1][j] - 2 * grid[i][j] + grid[i - 1][j]) / (dx ** 2);
            const d2Tdy2 = (grid[i][j + 1] - 2 * grid[i][j] + grid[i][j - 1]) / (dy ** 2);

            // 计算新温度并限制范围
            let newTemp = grid[i][j] + alpha * dt * (d2Tdx2 + d2Tdy2) + (Q[i][j] * dt);
            newTemp = Math.max(0, Math.min(150, newTemp));
            newGrid[i][j] = newTemp;
        }
    }

    // 替换原网格为新计算的网格
    grid = newGrid;
}

// 循环更新并渲染网格
function main() {
    updateGrid();
    renderGrid();
    requestAnimationFrame(main);
}

// 渲染网格到Canvas
function renderGrid() {
    ctx.clearRect(0, 0, canvas.width, canvas.height);

    for (let i = 0; i < width; i++) {
        for (let j = 0; j < height; j++) {
            // 约束hue范围在0-240之间
            let hue = Math.max(0, Math.min(240, ((100 - grid[i][j]) / 100) * 240));
            ctx.fillStyle = `hsl(${hue}, 100%, 50%)`;
            ctx.fillRect(i * dx, j * dy, dx, dy);
        }
    }
}

// 启动模拟
main();

内容的提问来源于stack exchange，提问作者Kirill Koshcheev