如何用Python求解该方程组并提取v₁、v₂、v₃、vₙ的数值？

我不清楚如何求解以下相互耦合的方程组，提取v₁、v₂、v₃和vₙ的数值用于打印计算。曾尝试用sympy符号求解，但无法将结果转换为int或float类型。

相关代码如下：

import cmath
import math
import numpy as np

r3 = math.sqrt(3)
pp = cmath.exp(-1j * math.pi * 2 / 3)
U = 400
U_ph = U / r3
U_a = U_ph
U_b = U_a * pp
U_c = U_b * pp

S_HP = 3e3 - cmath.sqrt((3e3 / 0.95) ** 2 - 3e3 ** 2) * 1j
S_PV = -3.6e3 + 0j
S_EV = 16 * 400 * r3 + 0j

Z = 0.3 + 0.1j
Z_HP = (U_ph**2 / S_HP).conjugate()
Z_PV = (U_ph**2 / S_PV).conjugate()
Z_EV = (U_ph**2 / S_EV).conjugate()

v_1 = ((U_a/Z) + (v_n/Z_HP)) / ((1/Z) + (1/Z_HP))
v_2 = ((U_b/Z) + (v_n/Z_PV)) / ((1/Z) + (1/Z_PV))
v_3 = ((U_c/Z) + (v_n/Z_EV)) / ((1/Z) + (1/Z_EV))
v_n = ((v_1/Z_HP) + (v_2/Z_PV) + (v_3/Z_EV)) / ((1/Z) + (1/Z_HP) + (1/Z_PV) + (1/Z_EV))

print(f"\n |v1-vn| = {np.abs(v_1 - v_n):.0f} V.")
print(f"\n |v2-vn| = {np.abs(v_2 - v_n):.0f} V.")
print(f"\n |v3-vn| = {np.abs(v_3 - v_n):.0f} V.")

解决方案 #

方法一：数值迭代法 #

由于v₁、v₂、v₃与vₙ相互依赖，可通过迭代逼近求解：

• 初始化v_n的初始值（比如0j）
• 循环计算v₁、v₂、v₃，再更新v_n，直到变量变化量小于设定阈值（比如1e-6）

代码示例：

import cmath
import math
import numpy as np

r3 = math.sqrt(3)
pp = cmath.exp(-1j * math.pi * 2 / 3)
U = 400
U_ph = U / r3
U_a = U_ph
U_b = U_a * pp
U_c = U_b * pp

S_HP = 3e3 - cmath.sqrt((3e3 / 0.95) ** 2 - 3e3 ** 2) * 1j
S_PV = -3.6e3 + 0j
S_EV = 16 * 400 * r3 + 0j

Z = 0.3 + 0.1j
Z_HP = (U_ph**2 / S_HP).conjugate()
Z_PV = (U_ph**2 / S_PV).conjugate()
Z_EV = (U_ph**2 / S_EV).conjugate()

# 迭代求解
v_n = 0j  # 初始值
threshold = 1e-6
max_iter = 1000
for _ in range(max_iter):
    # 计算当前v1, v2, v3
    v1_new = ((U_a/Z) + (v_n/Z_HP)) / ((1/Z) + (1/Z_HP))
    v2_new = ((U_b/Z) + (v_n/Z_PV)) / ((1/Z) + (1/Z_PV))
    v3_new = ((U_c/Z) + (v_n/Z_EV)) / ((1/Z) + (1/Z_EV))
    # 更新v_n
    v_n_new = ((v1_new/Z_HP) + (v2_new/Z_PV) + (v3_new/Z_EV)) / ((1/Z) + (1/Z_HP) + (1/Z_PV) + (1/Z_EV))
    # 检查收敛
    if np.abs(v_n_new - v_n) < threshold:
        v_n = v_n_new
        v_1 = v1_new
        v_2 = v2_new
        v_3 = v3_new
        break
    v_n = v_n_new

print(f"\n |v1-vn| = {np.abs(v_1 - v_n):.0f} V.")
print(f"\n |v2-vn| = {np.abs(v_2 - v_n):.0f} V.")
print(f"\n |v3-vn| = {np.abs(v_3 - v_n):.0f} V.")

方法二：SymPy符号求解并转换数值 #

用SymPy定义符号变量，联立方程求解后代入数值计算，注意要将符号解转换为复数数值：

代码示例：

import cmath
import math
import numpy as np
import sympy as sp

r3 = math.sqrt(3)
pp = cmath.exp(-1j * math.pi * 2 / 3)
U = 400
U_ph = U / r3
U_a = U_ph
U_b = U_a * pp
U_c = U_b * pp

S_HP = 3e3 - cmath.sqrt((3e3 / 0.95) ** 2 - 3e3 ** 2) * 1j
S_PV = -3.6e3 + 0j
S_EV = 16 * 400 * r3 + 0j

Z = 0.3 + 0.1j
Z_HP = (U_ph**2 / S_HP).conjugate()
Z_PV = (U_ph**2 / S_PV).conjugate()
Z_EV = (U_ph**2 / S_EV).conjugate()

# 定义符号变量
v1, v2, v3, vn = sp.symbols('v1 v2 v3 vn')

# 建立方程
eq1 = sp.Eq(v1, ((U_a/Z) + (vn/Z_HP)) / ((1/Z) + (1/Z_HP)))
eq2 = sp.Eq(v2, ((U_b/Z) + (vn/Z_PV)) / ((1/Z) + (1/Z_PV)))
eq3 = sp.Eq(v3, ((U_c/Z) + (vn/Z_EV)) / ((1/Z) + (1/Z_EV)))
eq4 = sp.Eq(vn, ((v1/Z_HP) + (v2/Z_PV) + (v3/Z_EV)) / ((1/Z) + (1/Z_HP) + (1/Z_PV) + (1/Z_EV)))

# 求解方程组
solution = sp.solve((eq1, eq2, eq3, eq4), (v1, v2, v3, vn))

# 提取数值（转换为复数）
v_1 = complex(solution[v1])
v_2 = complex(solution[v2])
v_3 = complex(solution[v3])
v_n = complex(solution[vn])

print(f"\n |v1-vn| = {np.abs(v_1 - v_n):.0f} V.")
print(f"\n |v2-vn| = {np.abs(v_2 - v_n):.0f} V.")
print(f"\n |v3-vn| = {np.abs(v_3 - v_n):.0f} V.")

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