书上例题的代码
Barczyk法
#Barczyk
from sympy import *
q1,r1 = 20,0.375
q2,r2 = 11.5,1.5
q3,r3 = 8.5,2
q4,r4 = 10.5,4.6875
q5,r5 = 1,12.5
q6,r6 = 9.5,1.7361
q7,r7 = 20,0.5
i = 0
while True:
i += 1
f1 = r5 * q5 * abs(q5) + r2 * q2 * abs(q2) - r3 * q3 * abs(q3)
f2 = -1 * r6 * q6 * abs(q6) + r2 * q2 * abs(q2) - r3 * q3 * abs(q3) + r4 * q4 * abs(q4)
df12 = 2*r5*abs(q5)+2*r2*abs(q2)+2*r3*abs(q3)
df14 = 0
df22 = 2*r6*abs(q6)+2*r2*abs(q2)+2*r3*abs(q3)+2*r4*abs(q4)
df24 = 2*r6*abs(q6)+2*r2*abs(q2)+2*r3*abs(q3)+2*r4*abs(q4)
print(f'第{i}次迭代:')
J = Matrix([[df12,df14],[df22,df24]])
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[0]))
max_f = max(abs(f1),abs(f2))
if max_f <= 10**-3 :
print(q1,q2,q3,q4,q5,q6)
print(f'△f = {max_f}')
break
else:
q2 = q2 + var_Q[0]
q4 = q4 + var_Q[1]
q3 = 20-q2
q5 = q2-q4
q6 = 20-q4
print(q1, q2, q3, q4,q5,q6)
print(f'△f = {max_f}')
结果:
第1次迭代:
J: Matrix([[93.5000000000000, 0], [199.923400000000, 199.923400000000]])
J的逆矩阵: Matrix([[0.0106951871657754, 0], [-0.0106951871657754, 0.00500191573372602]])
F: Matrix([[66.3750000000000], [413.988850000000]])
Q的变量 Matrix([[-0.709893048128342], [-1.36084429427380]])
20 10.7901069518717 9.20989304812834 9.13915570572620 1.65095124614545 10.8608442942738
△f = 413.98884999999996
第2次迭代:
J: Matrix([[110.483674201765, 0], [192.600501347889, 192.600501347889]])
J的逆矩阵: Matrix([[0.00905111101006476, 0], [-0.00905111101006476, 0.00519209448055240]])
F: Matrix([[39.0658523476833], [191.728107485435]])
Q的变量 Matrix([[-0.353589366301681], [-0.641881082340207]])
20 10.4365175855700 9.56348241443003 8.49727462338600 1.93924296218398 11.5027253766140
△f = 191.728107485435
第3次迭代:
J: Matrix([[118.044556469030, 0], [189.165195061353, 189.165195061353]])
J的逆矩阵: Matrix([[0.00847137750280220, 0], [-0.00847137750280219, 0.00528638473729623]])
F: Matrix([[27.4692480183935], [89.2076255559517]])
Q的变量 Matrix([[-0.232702369681913], [-0.238883460507508]])
20 10.2038152158881 9.79618478411194 8.25839116287849 1.94542405300957 11.7416088371215
△f = 89.2076255559517
第4次迭代:
J: Matrix([[118.431786109351, 0], [187.987816140351, 187.987816140351]])
J的逆矩阵: Matrix([[0.00844367912408813, 0], [-0.00844367912408813, 0.00531949368066175]])
F: Matrix([[11.5547291164030], [44.5905147077962]])
Q的变量 Matrix([[-0.0975644250246650], [-0.139634536180912]])
20 10.1062507908634 9.89374920913660 8.11875662669758 1.98749416416582 11.8812433733024
△f = 44.5905147077962
第5次迭代:
J: Matrix([[119.581103313282, 0], [187.261145825207, 187.261145825207]])
J的逆矩阵: Matrix([[0.00836252528445209, 0], [-0.00836252528445210, 0.00534013607357406]])
F: Matrix([[6.80857390257319], [21.3300428548339]])
Q的变量 Matrix([[-0.0569368714113290], [-0.0569684598886504]])
20 10.0493139194521 9.95068608054793 8.06178816680893 1.98752575264314 11.9382118331911
△f = 21.3300428548339
第6次迭代:
J: Matrix([[119.638829896626, 0], [186.981809271588, 186.981809271588]])
J的逆矩阵: Matrix([[0.00835849030673441, 0], [-0.00835849030673442, 0.00534811382933790]])
F: Matrix([[2.82899114806503], [10.6722292031993]])
Q的变量 Matrix([[-0.0236460950889390], [-