Python求解三元非线性方程组遇精度问题求助

我用Python求解含3个未知数（theta3、theta4、theta5）的3个非线性方程组时遇到了问题。尝试过差分进化法、basinhopping、minimize等多种优化方法，但求解得到的参数代入原方程后，计算出的Vcos2nu、Vsin2nu、Vdelta与给定值误差极大。

方程组 #

import math
from math import sin, cos, sqrt, pi, fabs

# 方程1：Vcos2nu表达式
Vcos2nu = (0.25*(-sin(2*theta4 + pi/4) + sqrt(2)*cos(2*theta4) - cos(2*theta4 + pi/4) + cos(-2*theta3 + 2*theta4 + pi/4) + cos(-2*theta4 + 2*theta5 + pi/4))**2 - 0.25*(-sin(2*theta4 + pi/4) + sqrt(2)*cos(2*theta4) - cos(2*theta4 + pi/4) + cos(2*theta3 - 2*theta4 + pi/4) + cos(2*theta4 - 2*theta5 + pi/4))**2)/(0.25*(-sin(2*theta4 + pi/4) + sqrt(2)*cos(2*theta4) - cos(2*theta4 + pi/4) + cos(-2*theta3 + 2*theta4 + pi/4) + cos(-2*theta4 + 2*theta5 + pi/4))**2 + 0.25*(-sin(2*theta4 + pi/4) + sqrt(2)*cos(2*theta4) - cos(2*theta4 + pi/4) + cos(2*theta3 - 2*theta4 + pi/4) + cos(2*theta4 - 2*theta5 + pi/4))**2)

# 方程2：Vsin2nu表达式
Vsin2nu = 2*fabs(-0.5*sin(2*theta4 + pi/4) + 0.5*sqrt(2)*cos(2*theta4) - 0.5*cos(2*theta4 + pi/4) + 0.5*cos(-2*theta3 + 2*theta4 + pi/4) + 0.5*cos(-2*theta4 + 2*theta5 + pi/4))*fabs(-0.5*sin(2*theta4 + pi/4) + 0.5*sqrt(2)*cos(2*theta4) - 0.5*cos(2*theta4 + pi/4) + 0.5*cos(2*theta3 - 2*theta4 + pi/4) + 0.5*cos(2*theta4 - 2*theta5 + pi/4))/(0.25*(-sin(2*theta4 + pi/4) + sqrt(2)*cos(2*theta4) - cos(2*theta4 + pi/4) + cos(-2*theta3 + 2*theta4 + pi/4) + cos(-2*theta4 + 2*theta5 + pi/4))**2 + 0.25*(-sin(2*theta4 + pi/4) + sqrt(2)*cos(2*theta4) - cos(2*theta4 + pi/4) + cos(2*theta3 - 2*theta4 + pi/4) + cos(2*theta4 - 2*theta5 + pi/4))**2)

# 方程3：Vdelta表达式（注意：原arg函数针对实数存在逻辑矛盾，需确认定义）
expr_delta = -cos(-2*theta3 + 2*theta4 + pi/4) + cos(2*theta3 - 2*theta4 + pi/4) - cos(-2*theta4 + 2*theta5 + pi/4) + cos(2*theta4 - 2*theta5 + pi/4)
Vdelta = math.atan2(0, expr_delta)  # 原arg(实数)仅能取0或π，与给定值不符，此处为临时写法

给定目标值 #

Vcos2nu = -0.28543
Vsin2nu = -0.13766479
Vdelta  = -0.94906616

问题分析与解决建议 #

1. 修正方程逻辑（核心问题） #

第三个方程的arg函数存在明显矛盾：数学中arg(z)是复数的辐角，若输入为实数，结果只能是0（正实数）或π（负实数），与给定的-0.949弧度完全不符。大概率是输入笔误，建议：

• 确认原方程是否应为Vdelta = expr_delta（直接让表达式结果等于目标值）
• 或是否是复数辐角的输入遗漏了虚部（比如arg(a + b*1j)）
• 或是否是atan2(expr_delta, 某个常数)的简写

2. 简化目标函数，降低计算误差 #

原表达式过于冗长，重复计算多，容易引入数值误差。先定义中间变量化简：

def objective(params):
    theta3, theta4, theta5 = params
    pi_val = math.pi
    # 定义中间变量减少重复计算
    A = -math.sin(2*theta4 + pi_val/4) + math.sqrt(2)*math.cos(2*theta4) - math.cos(2*theta4 + pi_val/4)
    B1 = math.cos(2*theta4 - 2*theta3 + pi_val/4)  # 等价于cos(-2θ3+2θ4+π/4)
    B2 = math.cos(2*theta5 - 2*theta4 + pi_val/4)  # 等价于cos(-2θ4+2θ5+π/4)
    C1 = math.cos(2*theta3 - 2*theta4 + pi_val/4)
    C2 = math.cos(2*theta4 - 2*theta5 + pi_val/4)
    
    sum1 = A + B1 + B2
    sum2 = A + C1 + C2
    denominator = sum1**2 + sum2**2
    if denominator < 1e-10:
        v_cos = 0.0
        v_sin = 0.0
    else:
        v_cos = (sum1**2 - sum2**2) / denominator
        v_sin = 2 * math.fabs(sum1) * math.fabs(sum2) / denominator
    
    # 假设Vdelta直接等于expr_delta（需根据实际方程修正）
    expr_delta = -B1 + C1 - B2 + C2
    
    # 计算误差平方和作为优化目标
    error = (v_cos + 0.28543)**2 + (v_sin + 0.13766479)**2 + (expr_delta + 0.94906616)**2
    return error

3. 强化全局搜索能力 #

非线性方程组易陷入局部最优，需调整优化策略：

• 使用differential_evolution时，扩大搜索范围（比如参数范围设为[(0, 2*math.pi), (0, 2*math.pi), (0, 2*math.pi)]），并调大popsize增强种群多样性
• 用basinhopping时，增加迭代次数与温度参数，提升跳出局部最优的概率
• 多设置不同初始点，多次运行后选择误差最小的结果

4. 验证代码示例 #

import math
from scipy.optimize import differential_evolution

# 目标值
TARGET_COS = -0.28543
TARGET_SIN = -0.13766479
TARGET_DELTA = -0.94906616

def objective(params):
    theta3, theta4, theta5 = params
    pi_val = math.pi
    A = -math.sin(2*theta4 + pi_val/4) + math.sqrt(2)*math.cos(2*theta4) - math.cos(2*theta4 + pi_val/4)
    B1 = math.cos(2*theta4 - 2*theta3 + pi_val/4)
    B2 = math.cos(2*theta5 - 2*theta4 + pi_val/4)
    C1 = math.cos(2*theta3 - 2*theta4 + pi_val/4)
    C2 = math.cos(2*theta4 - 2*theta5 + pi_val/4)
    
    sum1 = A + B1 + B2
    sum2 = A + C1 + C2
    denominator = sum1**2 + sum2**2
    v_cos = (sum1**2 - sum2**2)/denominator if denominator > 1e-10 else 0.0
    v_sin = 2 * math.fabs(sum1)*math.fabs(sum2)/denominator if denominator > 1e-10 else 0.0
    expr_delta = -B1 + C1 - B2 + C2
    
    return (v_cos-TARGET_COS)**2 + (v_sin-TARGET_SIN)**2 + (expr_delta-TARGET_DELTA)**2

# 参数范围（角度通常在0~2π之间）
bounds = [(0, 2*math.pi), (0, 2*math.pi), (0, 2*math.pi)]

# 差分进化求解
result = differential_evolution(objective, bounds, popsize=20, maxiter=1000)

# 输出结果
print(f"最优参数：theta3={result.x[0]:.6f}, theta4={result.x[1]:.6f}, theta5={result.x[2]:.6f}")
print(f"误差平方和：{result.fun:.6f}")

# 验证计算值
theta3, theta4, theta5 = result.x
A = -math.sin(2*theta4 + math.pi/4) + math.sqrt(2)*math.cos(2*theta4) - math.cos(2*theta4 + math.pi/4)
B1 = math.cos(2*theta4 - 2*theta3 + math.pi/4)
B2 = math.cos(2*theta5 - 2*theta4 + math.pi/4)
C1 = math.cos(2*theta3 - 2*theta4 + math.pi/4)
C2 = math.cos(2*theta4 - 2*theta5 + math.pi/4)
sum1 = A + B1 + B2
sum2 = A + C1 + C2
denominator = sum1**2 + sum2**2
v_cos = (sum1**2 - sum2**2)/denominator
v_sin = 2 * math.fabs(sum1)*math.fabs(sum2)/denominator
expr_delta = -B1 + C1 - B2 + C2
print(f"计算值：Vcos2nu={v_cos:.6f}, Vsin2nu={v_sin:.6f}, Vdelta={expr_delta:.6f}")
print(f"目标值：Vcos2nu={TARGET_COS:.6f}, Vsin2nu={TARGET_SIN:.6f}, Vdelta={TARGET_DELTA:.6f}")

内容的提问来源于stack exchange，提问作者H.AWAD