例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