流体网络拓扑（4）——网络分流算法中的Barczyk法和Cross法的例题

例1

代码：

#Barczyk，有回路附加阻力
from sympy import *

q1,r1 = 3,2
q2,r2 = 2,3
q3,r3 = 5,1


C_T = Matrix([[-1],[1]])
i = 0
while True:
    i += 1
    print(f'第{
     i}次迭代：')
    f1 = round(r1 * q1 * abs(q1) - r2 * q2 * abs(q2), 4)
    f2 = round(r3 * q3 * abs(q3) + r2 * q2 * abs(q2) - 10, 4)

    df11 = round(2*r1*abs(q1)+2*r2*abs(q2), 4)
    df13 = 0
    df21 = 0
    df23 = round(2*r3*abs(q3)+2*r2*abs(q2), 4)
    print('f1=',round(f1,4),'f2=',round(f2,4))
    print('df13=',df11,'df24=',df23)
    J = Matrix([[df11,df13],[df21,df23]])
    print('J:',J)
    J_inv = J**(-1)
    print('J的逆矩阵:',J_inv)
    F = Matrix([[f1],[f2]])
    print('F：',F)
    var_Q = -1*J_inv*F
    print('Q的变量',var_Q)
    max_vq = max(abs(var_Q[0]),abs(var_Q[1]))
    max_f = max(abs(f1),abs(f2))
    if max_vq <= 10**-6 and max_f <= 10**-6:
        print('q1=',q1, 'q2=',q2, 'q3=',q3)
        print(f'max△q = {
     round(max_vq,4)} , max△f = {
     round(max_f,4)}')
        break
    else:
        q1 = q1 + var_Q[0]
        q3 = q3 + var_Q[1]
        Q = Matrix([[q1,q3]]) * C_T
        q2 = Q[0]
        print('q1=',q1, 'q2=',q2, 'q3=',q3)
        print(f'max△q = {
     round(max_vq,4)} , max△f = {
     round(max_f,4)}')

结果：

第1次迭代：
f1= 6 f2= 27
df13= 24 df24= 22
J: Matrix([[24, 0], [0, 22]])
J的逆矩阵: Matrix([[1/24, 0], [0, 1/22]])
F： Matrix([[6], [27]])
Q的变量 Matrix([[-1/4], [-27/22]])
q1= 11/4 q2= 45/44 q3= 83/22
max△q = 1.2273 , max△f = 27
第2次迭代：
f1= 11.9871 f2= 7.3714
df13= 17.1364 df24= 13.6818
J: Matrix([[17.1364, 0], [0, 13.6818]])
J的逆矩阵: Matrix([[0.0583553, 0], [0, 0.0730898]])
F： Matrix([[11.9871], [7.3714]])
Q的变量 Matrix([[-0.699511], [-0.538774]])
q1= 2.05049 q2= 1.18346 q3= 3.23395
max△q = 0.6995 , max△f = 11.9871
第3次迭代：
f1= 4.2072 f2= 4.6602
df13= 15.3027 df24= 13.5687
J: Matrix([[15.3027, 0], [0, 13.5687]])
J的逆矩阵: Matrix([[0.0653479, 0], [0, 0.0736990]])
F： Matrix([[4.2072], [4.6602]])
Q的变量 Matrix([[-0.274932], [-0.343452]])
q1= 1.77556 q2= 1.11494 q3= 2.89050
max△q = 0.3435 , max△f = 4.6602
第4次迭代：
f1= 2.5759 f2= 2.0843
df13= 13.7919 df24= 12.4707
J: Matrix([[13.7919, 0], [0, 12.4707]])
J的逆矩阵: Matrix([[0.0725063, 0], [0, 0.0801880]])
F： Matrix([[2.5759], [2.0843]])
Q的变量 Matrix([[-0.186769], [-0.167136]])
q1= 1.58879 q2= 1.13458 q3= 2.72337
max△q = 0.1868 , max△f = 2.5759
第5次迭代：
f1= 1.1867 f2= 1.2785
df13= 13.1626 df24= 12.2542
J: Matrix([[13.1626, 0], [0, 12.2542]])
J的逆矩阵: Matrix([[0.0759728, 0], [0, 0.0816047]])
F： Matrix([[1.1867], [1.2785]])
Q的变量 Matrix([[-0.0901570], [-0.104332]])
q1= 1.49863 q2= 1.12040 q3= 2.61903
max△q = 0.1043 , max△f = 1.2785
第6次迭代：
f1= 0.7259 f2= 0.6252
df13= 12.7169 df24= 11.9605
J: Matrix([[12.7169, 0], [0, 11.9605]])
J的逆矩阵: Matrix([[0.0786355, 0], [0, 0.0836085]])
F： Matrix([[0.7259], [0.6252]])
Q的变量 Matrix([[-0.0570814], [-0.0522719]])
q1= 1.44155 q2= 1.12521 q3= 2.56676
max△q = 0.0571 , max△f = 0.7259
第7次迭代：
f1= 0.3578 f2= 0.3866
df13= 12.5175 df24= 11.8848
J: Matrix([[12.5175, 0], [0, 11.8848]])
J的逆矩阵: Matrix([[0.0798882, 0], [0, 0.0841411]])
F： Matrix([[0.3578], [0.3866]])
Q的变量 Matrix([[-0.0285839], [-0.0325290]])
q1= 1.41297 q2= 1.12127 q3= 2.53423
max△q = 0.0325 , max△f = 0.3866
第8次迭代：
f1= 0.2212 f2= 0.1941
df13= 12.3795 df24= 11.7961
J: Matrix([[12.3795, 0], [0, 11.7961]])
J的逆矩阵: Matrix([[0.0807787, 0], [0, 0.0847738]])
F： Matrix([[0.2212], [0.1941]])
Q的变量 Matrix([[-0.0178683], [-0.0164545]])
q1= 1.39510 q2= 1.12268 q3= 2.51778
max△q = 0.0179 , max△f = 0.2212
第9次迭代：
f1= 0.1114 f2= 0.1204
df13= 12.3165 df24= 11.7716
J: Matrix([[12.3165, 0], [0, 11.7716]])
J的逆矩阵: Matrix([[0.0811919, 0], [0, 0.0849502]])
F： Matrix([[0.1114], [0.1204]])
Q的变量 Matrix([[-0.00904475], [-0.0102280]])
q1= 1.38605 q2= 1.12150 q3= 2.50755
max△q = 0.0102 , max△f = 0.1204
第10次迭代：
f1= 0.0690 f2=