梯度法三维Python代码实现这是⼀次反演课作业，新⼿编写过程有点繁琐。
"""
最速下降法中间穿插割线法
"""
from numpy import*
from sympy import symbols, diff
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
import matplotlib as mpl
x1, x2, x3, x4 = symbols('x1, x2, x3, x4')
x_0 =[0,0]
print('x0=', x_0)
alf =[x4]
f =100*(x2-x1**2)**2+(1-x1)**2
k1 =[]
k2 =[]
k3 =[]
F =[]
F_X0 =[]
df =[diff(f, x1), diff(f, x2)]# 求⼀阶导
'''
求df(x0) 进⼊迭代
'''
def steepset(xx1, xx2):
x0 =[xx1, xx2]
# print('x0=', x0)
g_0 =[]
for i in range(0,2):# i 的取值范围与⾃变量有关下同
aa = df[i]
# g_0.append(aa.subs({x1: 4, x2: 2, x3: -1}))
g_0.append(aa.subs({x1: xx1, x2: xx2}))# 不同⾃变量输⼊不同
print('g=', g_0)
a2 =[]
a1 = mat(x0)- mat(alf)*mat(g_0)# <class 'numpy.matrixlib.defmatrix.matrix'>
a1 = np.array(a1)# 将复杂矩阵转化为简单矩阵  <class 'numpy.ndarray'>
print('a1=', a1)
for i in range(0,len(a1)):
for j in range(0,2):
a2.append(a1[i][j])
print('a2=', a2)
sef = f.subs({x1: a2[0], x2: a2[1]})
print('sef=', sef)
"""
割线法，求α
"""
# x = symbols('x')
# g = x**3 - 12.2*x**2 + 7.45*x + 42    # 低次使⽤ g
# g = x**3 - 12.2*x**2 + 7.45*x + 42    # 低次使⽤ g
# y = [13, 11]
# sef = 4*(4-1024*x)**4 + (2*x-1)**2        # ⾼次使⽤导数
y =[0.02,0.0001]# 割线法起始点，y只是⼀个代号起点很关键，影响后续计算
#
def func(a, b):
# 使⽤导数形式
sex_k1 =(diff(sef).subs({x4: a})*b - diff(sef).subs({x4: b})*a)/(diff(sef).subs({x4: a})- diff(sef).subs({x4: b})) return sex_k1
# 使⽤g函数
# x_k1 = (g.subs({x: a}) * b - g.subs({x: b}) * a) / (g.subs({x: a}) - g.subs({x: b}))
# return x_k1
# y.append(func(y[0], y[1]))
# print(y)
#
#
n =0
for i in range(0,7):# 迭代次数⾃⼰选择
y.append(func(y[i], y[i+1]))
n +=1
# print('第{}迭代的结果为：{}'.format(n, func(y[i], y[i+1])))
alf_1 =[]
for i in range(0,len(y)):
if str(y[-1])=='nan':
if str(y[i])=='nan':
python新手代码画图alf_1.append(y[i -1])
# print('经过{}次迭代的结果为：{}'.format(n - 2, alf_1[0]))
break
# break
if str(y[-1])!='nan':
alf_1.append(y[-1])
# print('经过{}次迭代的结果为：{}'.format(n, alf_1[0]))
break
print('alf=', alf_1[0])
xn =[]
x = mat(x0).T - alf_1[0]*mat(g_0).T
x = np.array(x)
for i in range(0,len(x)):
xn.append(x[i][0])
# print('x=', x[i][0])
# print(type(x[i][0]))
return xn
'''
以上求出α，下⼀步求x1，进⽽求df(x1), 使⽤迭代求df(xn)
'''
x_3 =[x_0[0]]
x_4 =[x_0[1]]
m =0
while m <200:# 迭代次数⾃⼰选择
x_0 = steepset(x_0[0], x_0[1])
m +=1
if str(x_0[0])=='zoo':
break
x_3.append(x_0[0])
x_4.append(x_0[1])
print("x{}={}\n".format(m,  x_0))
print('f=', f.subs({x1: x_0[0], x2: x_0[1]}))
print('x_0=', x_0)
print('x_0=', x_0)
print('x_3=', x_3)
print('x_4=', x_4)
#下两⾏引⽤⾃ blog.csdn/yiyue21/article/details/89574413?utm_source=Params['font.sans-serif']=[u'SimHei']# 解决标题中⽂⽆法显⽰问题
fig = plt.figure()
ax = Axes3D(fig)
x1 = np.arange(-1,1.5,0.05)
x2 = np.arange(-1,1.5,0.05)
def function(x3, x4):
return100*(x4-x3**2)**2+(1-x3)**2
X, Y = np.meshgrid(x1, x2)# ⽹格的创建，这个是关键
Z =100*(Y-X**2)**2+(1-X)**2
print('x_3=', x_3)
print('x_4=', x_4)
x3 = np.array(x_3)
x4 = np.array(x_4)
x1_list =[]
x2_list =[]
y_list =[]
for i in range(0,len(x_3)):# 将梯度所⽤的点收集⽤于画图
x1_list.append(x_3[i])
x2_list.append(x_4[i])
y_list.append(function(x_3[i], x_4[i]))
print('f=', y_list)
plt.xlabel('x1')
plt.ylabel('x2')
ax.plot(x1_list, x2_list, y_list,'r*--')# 画出梯度线
ax.set_title('起点(0,0)割线法起点(0.02, 0.0001)迭代33次 f=5e-25')# 设置标题
plt.show()
运⾏结果如下

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。