Python函数返回非物理值问题排查：ηₖ计算异常求助

函数定义 #

给定递归数学函数：

\eta_{k}=\frac{\epsilon*2^{k-1}}{\sin\theta} \cdot \left(1-\cos\theta \prod_{j=1}^{k-1}(1+\sqrt{1-\eta_{j}})\right)
\quad \text{for } k>1
\\
\eta_{1}=\epsilon \cdot \tan\left(\frac{\theta}{2}\right)

其中：

• k 为整数
• \epsilon 和 \theta 为浮点数
• 根据定义，\eta_{k} 的取值不应超过 2^{k-1}

问题现象 #

调用函数 eta(1.0000001, 3, np.pi/10000000) 时，得到结果：

Decimal('851903.0534777333315363710505')

该结果远大于 2^{3-1}=4，属于不符合定义的非物理值。

待排查代码 #

from decimal import *
from math import *
import matplotlib.pyplot as plt
import numpy as np

getcontext().prec = 28


def eta(epsilon, k, theta):
    if k == 1:
        numerator = Decimal(epsilon) * Decimal((Decimal(1) - Decimal(cos(theta))))
        denominator = Decimal(sin(theta))
        return Decimal(numerator / denominator)
    else:
        product_term = 1
        for j in range(1, k):
            eta_j = Decimal(eta(epsilon, j, theta))
            product_term *= Decimal((Decimal(1) + Decimal(sqrt(1 - eta_j))))
        numerator = Decimal(epsilon)*Decimal(2**(k-1))*(Decimal(1)-Decimal(cos (theta))*Decimal(product_term)/Decimal(2**(k-1)))
        denominator = Decimal(sin(theta)) 
        return numerator / denominator

错误分析与修正 #

1. 核心公式错误 #

代码中错误地在公式的括号内给 cosθ * product_term 添加了除以 2^{k-1} 的操作，而原公式中并无此步骤。这是导致结果严重偏离预期的根本原因：

• 原公式括号内为：1 - cosθ \cdot \prod_{j=1}^{k-1}(1+\sqrt{1-\eta_{j}})
• 代码中写成了：1 - \frac{cosθ \cdot product_term}{2^{k-1}}

2. 精度相关问题 #

• 使用math.sqrt处理Decimal类型数值时，会先将Decimal转为float，丢失精度，应改用Decimal类型自带的.sqrt()方法。
• Decimal(2**(k-1)) 先计算整数幂再转Decimal，当k较大时可能溢出或丢失精度，应改为Decimal(2) ** (k-1)直接进行Decimal幂运算。

修正后的代码 #

from decimal import *
import matplotlib.pyplot as plt
import numpy as np

getcontext().prec = 28


def eta(epsilon, k, theta):
    eps = Decimal(epsilon)
    theta_dec = Decimal(theta)
    cos_theta = theta_dec.cos()
    sin_theta = theta_dec.sin()
    
    if k == 1:
        # tan(theta/2) = (1 - cosθ)/sinθ，与原定义等价
        tan_half = (Decimal(1) - cos_theta) / sin_theta
        return eps * tan_half
    else:
        product_term = Decimal(1)
        for j in range(1, k):
            eta_j = eta(epsilon, j, theta)
            sqrt_term = (Decimal(1) - eta_j).sqrt()
            product_term *= (Decimal(1) + sqrt_term)
        # 严格按照原公式计算
        power_of_two = Decimal(2) ** (k-1)
        numerator = eps * power_of_two * (Decimal(1) - cos_theta * product_term)
        return numerator / sin_theta

修正后调用原测试用例，结果会符合η_k ≤ 2^{k-1}的定义要求。

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