Python求解给定x、y范围下非线性方程的z值

首先把你的方程整理为关于z的一元二次方程，用求根公式直接计算比数值求解器更高效。

方程整理 #

原方程：
z*cos(theta) - x*sin(theta) = A - B*[z*sin(theta) + x*cos(theta)]² - B*y²

令C = cos(theta)，S = sin(theta)，展开并整理成a*z² + b*z + c = 0的形式：

• 二次项系数：a = B * S²
• 一次项系数：b = C + 2*B*S*C*x
• 常数项：c = -x*S - A + B*x²*C² + B*y²

当判别式delta = b² - 4*a*c ≥ 0时，存在实根z：
z = [-b ± sqrt(delta)] / (2*a)

Python实现代码 #

import numpy as np

# 替换为你的实际参数值
theta = np.pi/4  # 示例角度，可修改
A = 10
B = 0.01

# 生成x和y的网格，步长可按需调整
x = np.linspace(-40, 40, 81)  # 从-40到40共81个点
y = np.linspace(-40, 40, 81)
X, Y = np.meshgrid(x, y)

# 预计算cos和sin值
C = np.cos(theta)
S = np.sin(theta)

# 计算二次方程的三个系数
a = B * S**2
b = C + 2 * B * S * C * X
c = -X * S - A + B * (X**2 * C**2 + Y**2)

# 计算判别式，筛选有实根的区域
delta = b**2 - 4 * a * c
valid_mask = delta >= 0

# 初始化z的两个根数组，无实根的地方设为NaN
z1 = np.full_like(X, np.nan)
z2 = np.full_like(X, np.nan)

# 计算有实根处的z值
sqrt_delta = np.sqrt(delta[valid_mask])
z1[valid_mask] = (-b[valid_mask] + sqrt_delta) / (2 * a)
z2[valid_mask] = (-b[valid_mask] - sqrt_delta) / (2 * a)

# 输出结果示例：比如取x=0,y=0处的z值
idx_x = np.where(x == 0)[0][0]
idx_y = np.where(y == 0)[0][0]
print(f"x=0,y=0时，z的实根为：{z1[idx_y, idx_x]}, {z2[idx_y, idx_x]}")

# 若需要保存所有z值到文件
np.savez('z_values.npz', x=x, y=y, z1=z1, z2=z2)

补充说明 #

• 如果a=0（即sin(theta)=0，theta为0或π），方程退化为一元一次方程，单独处理：z = (A - B*x² - B*y² + x*S) / C（此时S=0，C=±1）
• 代码中用NaN标记无实根的位置，方便后续筛选或可视化
• 可调整linspace的步长控制计算精度，步长越小结果越精细

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