RuntimeWarning: overflow encountered in double

RuntimeWarning: overflow encountered in double_scalars问题咨询

阿华AIGC实验室

2026-5-7

解决RuntimeWarning: overflow encountered in double_scalars与inf计算结果问题

首先帮你拆解这个警告：RuntimeWarning: overflow encountered in double_scalars 指的是双精度浮点数的标量运算中（Python里默认的float就是双精度类型，double_scalars就是单个双精度浮点数），计算出的数值超出了双精度能表示的最大范围（大概是±1.797e+308），直接溢出成了inf（无穷大）。

接下来看你的代码，核心问题其实有两个：

1. 列表引用赋值导致的迭代错误

你代码里的H_new = H_old是个大坑！这不是创建新的列表，只是把H_new指向了和H_old同一个内存对象。也就是说，你在循环里修改H_new[xind]的时候，H_old的元素也会跟着变——后面计算用的H_old已经是被更新过的值，不是上一个时间步的原始状态了。这种错误会让数值在迭代里疯狂累积，很快就超出双精度范围引发溢出。

2. 数值计算的稳定性隐患

另外，你用普通Python列表处理数值运算，不仅效率低，还容易出现引用问题；加上你的通量计算（v1、v2）涉及H**4这种高次幂，一旦数值因为迭代错误异常增长，分分钟就溢出。

针对性解决方案

先修复最致命的引用问题：改用numpy数组

把H_old和H_new换成numpy数组，用copy()创建副本，彻底切断引用关联：

# 初始化部分直接用numpy广播创建数组，不用手动循环
xint = np.arange(-L, L+dx, dx)
H_old = H_0 * (1 - np.abs(xint)/L)
H_new = H_old.copy()  # 明确创建副本，不是引用
result = [H_old.copy()]  # 存副本，避免后续修改影响历史结果

修正循环边界，避免索引越界

原代码里range(1, len(xint))会遍历到最后一个索引，虽然x=L时会进入else分支，但为了保险，把循环范围改成range(1, len(xint)-1)，确保xind+1不会超出数组范围，再单独处理首尾边界：

for xind in range(1, len(xint)-1):
    x = xint[xind]
    if abs(x) < L:  # 简化判断逻辑
        dH1 = (H_old[xind]-H_old[xind-1])/dx
        H1 = np.clip(H_old[xind], 0, H_0)  # 限制H的范围，避免负数/过大值
        dH2 = (H_old[xind+1]-H_old[xind])/dx
        H2 = np.clip(H_old[xind+1], 0, H_0)
        v1 = - (2/4)/(n+1) * fact * ((abs(dH1))**(n-1)) * dH1 * H1**4
        v2 = - (2/4)/(n+1) * fact * ((abs(dH2))**(n-1)) * dH2 * H2**4
        dH_dt = - (H2*v2-H1*v1)/dx + a
        H_new[xind] = H_old[xind] + dt * dH_dt
    else:
        H_new[xind] = 0.0
# 单独处理首尾两个边界点
H_new[0] = 0.0
H_new[-1] = 0.0

优化数值稳定性

如果还是溢出，可以试试减小时间步长dt（比如改成0.1），过大的dt会让数值迭代不稳定，容易炸掉。
用np.clip()限制H的范围，避免H**4因为数值异常变大直接溢出。

修改后的完整代码

## IMPORT PACKAGES
import numpy as np
import matplotlib.pyplot as plt

## DEFINITION AND VALUE OF PARAMETERS
A = 10**-16
rho = 910
g = 9.81
a = 1
n = 3
fact = A * ((rho * g)**n)
t_end = 10
dt = 1
dx = 100
L = 300000
b = 0
H_0 = 1000

## Solve equation numerically
# initialization
xint = np.arange(-L, L+dx, dx)
H_old = H_0 * (1 - np.abs(xint)/L)
H_new = H_old.copy()
result = [H_old.copy()]
print(len(xint), len(H_old))

# Loops through the time-span
for t in range(0, t_end+dt, dt):
    print(t, np.min(H_old), np.max(H_old))
    # Loops through the x-range
    for xind in range(1, len(xint)-1):
        x = xint[xind]
        if abs(x) < L:
            dH1 = (H_old[xind]-H_old[xind-1])/dx
            H1 = np.clip(H_old[xind], 0, H_0)
            dH2 = (H_old[xind+1]-H_old[xind])/dx
            H2 = np.clip(H_old[xind+1], 0, H_0)
            v1 = - (2/4)/(n+1) * fact * ((abs(dH1))**(n-1)) * dH1 * H1**4
            v2 = - (2/4)/(n+1) * fact * ((abs(dH2))**(n-1)) * dH2 * H2**4
            dH_dt = - (H2*v2-H1*v1)/dx + a
            H_new[xind] = H_old[xind] + dt * dH_dt
        else:
            H_new[xind] = 0.0
    # 处理边界
    H_new[0] = 0.0
    H_new[-1] = 0.0
    # 更新旧数组并存储结果
    H_old = H_new.copy()
    result.append(H_old.copy())