问题:
给定若干个三维离散点,任意求出一个方程使得尽可能的将这些点包含在其中,并用python画出函数图像
思路:
最小二乘法,根据最优性表示出一个表达式E,为了让E最小,那么对任意一个系数求偏导后的值应该为0,列出若干条式子后,再用内置的高斯消元,求解系数。百度文库中有若干篇论文描述最小二乘法的实现方法。
代码:
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
def fun(x): #处理符号问题
round(x,2)
if x>=0: return '+'+str(x)
else: return str(x)
#主函数
if __name__ == '__main__':
n=int(input("请输入一个n,代表点的个数"))
#print(n)
X=[]
Y=[]
Z=[]
#读入n个点的坐标
for i in range(n):
a=[int(i) for i in input('please input: ').split()]
X.append(a[0])
Y.append(a[1])
Z.append(a[2])
#print(X)
#print(Y)
#print(Z)
#fig = plt.figure() #建立一个新的图像
#ax = Axes3D(fig) #建立一个三维坐标系
#ax.scatter(X,Y,Z,c='b')
#ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
# plt.show()
#求方程系数
sigma_x = 0
for i in X : sigma_x+=i
sigma_y = 0
for i in Y : sigma_y+=i
sigma_z = 0
for i in Z : sigma_z+=i
sigma_x2 = 0
for i in X : sigma_x2+=i*i
sigma_y2 = 0
for i in Y : sigma_y2+=i*i
sigma_x3 = 0
for i in X : sigma_x3+=i*i*i
sigma_y3 = 0
for i in Y : sigma_y3+=i*i*i
sigma_x4 = 0
for i in X : sigma_x4+=i*i*i*i
sigma_y4 = 0
for i in Y : sigma_y4+=i*i*i*i
sigma_x_y = 0
for i in range(n) :
sigma_x_y+=X[i]*Y[i]
#print(sigma_xy)
sigma_x_y2 = 0
for i in range(n) : sigma_x_y2+=X[i]*Y[i]*Y[i]
sigma_x_y3 = 0
for i in range(n) : sigma_x_y3+=X[i]*Y[i]*Y[i]*Y[i]
sigma_x2_y = 0
for i in range(n) : sigma_x2_y+=X[i]*X[i]*Y[i]
sigma_x2_y2 = 0
for i in range(n) : sigma_x2_y2+=X[i]*X[i]*Y[i]*Y[i]
sigma_x3_y = 0
for i in range(n) : sigma_x3_y+=X[i]*X[i]*X[i]*Y[i]
sigma_z_x2 = 0
for i in range(n) : sigma_z_x2+=Z[i]*X[i]*X[i]
sigma_z_y2 = 0
for i in range(n) : sigma_z_y2+=Z[i]*Y[i]*Y[i]
sigma_z_x_y = 0
for i in range(n) : sigma_z_x_y+=Z[i]*X[i]*Y[i]
sigma_z_x = 0
for i in range(n) : sigma_z_x+=Z[i]*X[i]
sigma_z_y = 0
for i in range(n) : sigma_z_y+=Z[i]*Y[i]
#print("-----------------------")
#给出对应方程的矩阵形式
a=np.array([[sigma_x4,sigma_x3_y,sigma_x2_y2,sigma_x3,sigma_x2_y,sigma_x2],
[sigma_x3_y,sigma_x2_y2,sigma_x_y3,sigma_x2_y,sigma_x_y2,sigma_x_y],
[sigma_x2_y2,sigma_x_y3,sigma_y4,sigma_x_y2,sigma_y3,sigma_y2],
[sigma_x3,sigma_x2_y,sigma_x_y2,sigma_x2,sigma_x_y,sigma_x],
[sigma_x2_y,sigma_x_y2,sigma_y3,sigma_x_y,sigma_y2,sigma_y],
[sigma_x2,sigma_x_y,sigma_y2,sigma_x,sigma_y,n]])
b=np.array([sigma_z_x2,sigma_z_x_y,sigma_z_y2,sigma_z_x,sigma_z_y,sigma_z])
#高斯消元解线性方程
res= np.linalg.solve(a,b)
#print(a)
#print(b)
#print(x)
# print("-----------------------")
#输出方程形式
print("z=%.6s*x^2%.6s*xy%.6s*y^2%.6s*x%.6s*y%.6s"%(fun(res[0]),fun(res[1]),fun(res[2]),fun(res[3]),fun(res[4]),fun(res[5])))
#画曲面图和离散点
fig = plt.figure()#建立一个空间
ax =fig.add_subplot(111,projection='3d')# 3D坐标
n = 256
u = np.linspace(-20,20,n) # 创建一个等差数列
x,y = np.meshgrid(u,u) #转化成矩阵
#给出方程
z=res[0]*x*x+res[1]*x*y+res[2]*y*y+res[3]*x+res[4]*y+res[5]
#画出曲面
ax.plot_surface(x,y,z,rstride=3,cstride=3,cmap=cm.jet)
#画出点
ax.scatter(X, Y, Z,c='r')
参考网站:
https://blog.csdn.net/qq_39735236/article/details/79066032
https://www.cnblogs.com/NanShan2016/p/5493429.html
http://www.cnblogs.com/jukan/p/7826851.html