如下图所示，第一行是X，第一列是Y，其余为Z，欲拟合Z（X，Y）的二次关系式。
matlab代码
注意，要将X 和Y 扩充至与Z同维！！！！
%matlab
clc;clear;
x=[0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000. 0. 1000. 1285.714286 1571.428571 1857.142857 2142.857143 2428.571429 2714.285714 3000. 3285.714286 3571.428571 3857.142857 4142.857143 4428.571429 4714.285714 5000.]';
y=[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 40.17 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 212.6 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 373.7 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 522.2 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 658.8 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 787.7 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 914.4 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1042. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1168. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1293. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1414. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1534. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1657. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1784. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 1917. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2056. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2197. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2335. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467. 2467.]';
z=[0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.968730e-08 1.114980e-07 2.088450e-07 2.821970e-07 3.367890e-07 3.859630e-07 4.437040e-07 5.163040e-07 6.189590e-07 7.660220e-07 0.000000e+00 8.657850e-05 1.120380e-04 1.374110e-04 1.627070e-04 1.879460e-04 2.131280e-04 2.382670e-04 2.633660e-04 2.884420e-04 3.135050e-04 3.385640e-04 3.636300e-04 3.887080e-04 4.138110e-04 4.389530e-04 0.000000e+00 1.747260e-04 2.251530e-04 2.755080e-04 3.258000e-04 3.760330e-04 4.262250e-04 4.763780e-04 5.265080e-04 5.766110e-04 6.267010e-04 6.767880e-04 7.268780e-04 7.769820e-04 8.271050e-04 8.772580e-04 0.000000e+00 2.627510e-04 3.381570e-04 4.135090e-04 4.888110e-04 5.640640e-04 6.392810e-04 7.144710e-04 7.896350e-04 8.647800e-04 9.399170e-04 1.015052e-03 1.090189e-03 1.165335e-03 1.240493e-03 1.315672e-03 0.000000e+00 3.506750e-04 4.510800e-04 5.514370e-04 6.517510e-04 7.520420e-04 8.523030e-04 9.525350e-04 1.052752e-03 1.152959e-03 1.253158e-03 1.353349e-03 1.453540e-03 1.553742e-03 1.653950e-03 1.754176e-03 0.000000e+00 4.385230e-04 5.639250e-04 6.892960e-04 8.146460e-04 9.399770e-04 1.065288e-03 1.190586e-03 1.315869e-03 1.441141e-03 1.566409e-03 1.691672e-03 1.816933e-03 1.942196e-03 2.067466e-03 2.192754e-03 0.000000e+00 5.263150e-04 6.767290e-04 8.271260e-04 9.775100e-04 1.127889e-03 1.278260e-03 1.428622e-03 1.578978e-03 1.729328e-03 1.879673e-03 2.030012e-03 2.180346e-03 2.330685e-03 2.481029e-03 2.631377e-03 0.000000e+00 6.140900e-04 7.895090e-04 9.649330e-04 1.140361e-03 1.315790e-03 1.491221e-03 1.666652e-03 1.842083e-03 2.017513e-03 2.192939e-03 2.368362e-03 2.543784e-03 2.719205e-03 2.894625e-03 3.070050e-03 0.000000e+00 7.018540e-04 9.022860e-04 1.102738e-03 1.303204e-03 1.503686e-03 1.704186e-03 1.904689e-03 2.105197e-03 2.305706e-03 2.506216e-03 2.706724e-03 2.907227e-03 3.107729e-03 3.308231e-03 3.508735e-03 0.000000e+00 7.896490e-04 1.015085e-03 1.240558e-03 1.466066e-03 1.691599e-03 1.917156e-03 2.142726e-03 2.368309e-03 2.593901e-03 2.819490e-03 3.045081e-03 3.270675e-03 3.496263e-03 3.721850e-03 3.947432e-03 0.000000e+00 8.774890e-04 1.127933e-03 1.378421e-03 1.628958e-03 1.879538e-03 2.130144e-03 2.380778e-03 2.631431e-03 2.882094e-03 3.132766e-03 3.383442e-03 3.634115e-03 3.884785e-03 4.135453e-03 4.386112e-03 0.000000e+00 9.654110e-04 1.240842e-03 1.516335e-03 1.791894e-03 2.067506e-03 2.343163e-03 2.618853e-03 2.894566e-03 3.170297e-03 3.446041e-03 3.721793e-03 3.997543e-03 4.273291e-03 4.549033e-03 4.824764e-03 0.000000e+00 1.053412e-03 1.353825e-03 1.654321e-03 1.954894e-03 2.255529e-03 2.556214e-03 2.856947e-03 3.157717e-03 3.458507e-03 3.759315e-03 4.060135e-03 4.360956e-03 4.661773e-03 4.962585e-03 5.263383e-03 0.000000e+00 1.141520e-03 1.466917e-03 1.792393e-03 2.117965e-03 2.443609e-03 2.769320e-03 3.095085e-03 3.420891e-03 3.746728e-03 4.072588e-03 4.398463e-03 4.724344e-03 5.050224e-03 5.376096e-03 5.701943e-03 0.000000e+00 1.229768e-03 1.580117e-03 1.930574e-03 2.281130e-03 2.631766e-03 2.982480e-03 3.333262e-03 3.684096e-03 4.034967e-03 4.385866e-03 4.736786e-03 5.087715e-03 5.438638e-03 5.789543e-03 6.140437e-03 0.000000e+00 1.318161e-03 1.693457e-03 2.068870e-03 2.444395e-03 2.820015e-03 3.195721e-03 3.571499e-03 3.947342e-03 4.323233e-03 4.699154e-03 5.075095e-03 5.451048e-03 5.827001e-03 6.202938e-03 6.578851e-03 0.000000e+00 1.406731e-03 1.806947e-03 2.207302e-03 2.607777e-03 3.008360e-03 3.409038e-03 3.809803e-03 4.210636e-03 4.611521e-03 5.012443e-03 5.413394e-03 5.814359e-03 6.215315e-03 6.616257e-03 7.017177e-03 0.000000e+00 1.495475e-03 1.920609e-03 2.345884e-03 2.771293e-03 3.196820e-03 3.622454e-03 4.048181e-03 4.473980e-03 4.899840e-03 5.325747e-03 5.751682e-03 6.177630e-03 6.603573e-03 7.029502e-03 7.455402e-03 0.000000e+00 1.584421e-03 2.034449e-03 2.484634e-03 2.934958e-03 3.385412e-03 3.835981e-03 4.286642e-03 4.737389e-03 5.188202e-03 5.639061e-03 6.089956e-03 6.540868e-03 6.991775e-03 7.442662e-03 7.893517e-03 0.000000e+00 1.673583e-03 2.148484e-03 2.623560e-03 3.098776e-03 3.574132e-03 4.049610e-03 4.525192e-03 5.000857e-03 5.476601e-03 5.952399e-03 6.428226e-03 6.904072e-03 7.379918e-03 7.855743e-03 8.331530e-03]';
c=ones(length(x),1);
[X,Y]=meshgrid(linspace(0,50,30),linspace(0,50,30));
C=ones(30);
b=regress(z,[x.^2,y.^2,x.*y,x,y,c]);%添加非线性项进行拟合
figure,plot3(x,y,z,'p');
hold on;
mesh(X,Y,b(1)*X.^2+b(2)*Y.^2+b(3)*X.*Y+b(4)*X+b(5)*Y+b(6)*C);
disp(b(1));
disp(b(2));
disp(b(3));
disp(b(4));
disp(b(5));
disp(b(6));
结果：
-4.4377e-14 2.5213e-11 6.7171e-10 1.4773e-08 -5.9435e-08 2.1294e-05
------------------------------------分割线-----------------------------------
python代码
# python
from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import csv
def read_csv():
fuel_data_path = "data.csv"
with open(fuel_data_path, 'r') as f:
reader = csv.reader(f)
fuel_map = np.array(list(reader), dtype=float) #第一行是X，第一列是Y
# print(fuel_map)
X = fuel_map[1:2,1:]
X= np.tile(X, 20)
X= np.reshape(X,(20,16))
Y = fuel_map[2:,0]
Y = np.repeat(Y, 16)
Y = np.reshape(Y, (20, 16))
Z = fuel_map[2:,1:]
# print("X",np.shape(X))
# print("Y",np.shape(Y))
# print("Z",np.shape(Z))
return X,Y,Z
def fun(x):
round(x, 2)
if x >= 0:
return '+'+str(x)
else:
return str(x)
def get_res(X, Y, Z, n):
# 求方程系数
sigma_x = 0
for i in range(n): sigma_x += X[0, i]
sigma_y = 0
for i in range(n): sigma_y += Y[0, i]
sigma_z = 0
for i in range(n): sigma_z += Z[0, i]
sigma_x2 = 0
for i in range(n): sigma_x2 += X[0, i] * X[0, i]
sigma_y2 = 0
for i in range(n): sigma_y2 += Y[0, i] * Y[0, i]
sigma_x3 = 0
for i in range(n): sigma_x3 += X[0, i] * X[0, i] * X[0, i]
sigma_y3 = 0
for i in range(n): sigma_y3 += Y[0, i] * Y[0, i] * Y[0, i]
sigma_x4 = 0
for i in range(n): sigma_x4 += X[0, i] * X[0, i] * X[0, i] * X[0, i]
sigma_y4 = 0
for i in range(n): sigma_y4 += Y[0, i] * Y[0, i] * Y[0, i] * Y[0, i]
sigma_x_y = 0
for i in range(n): sigma_x_y += float(X[0, i] * Y[0, i])
sigma_x_y2 = 0
for i in range(n): sigma_x_y2 += float(X[0, i] * Y[0, i] * Y[0, i])
sigma_x_y3 = 0
for i in range(n): sigma_x_y3 += float(X[0, i] * Y[0, i] * Y[0, i] * Y[0, i])
sigma_x2_y = 0
for i in range(n): sigma_x2_y += float(X[0, i] * X[0, i] * Y[0, i])
sigma_x2_y2 = 0
for i in range(n): sigma_x2_y2 += float(X[0, i] * X[0, i] * Y[0, i] * Y[0, i])
sigma_x3_y = 0
for i in range(n): sigma_x3_y += float(X[0, i] * X[0, i] * X[0, i] * Y[0, i])
sigma_z_x2 = 0
for i in range(n): sigma_z_x2 += float(Z[0, i] * X[0, i] * X[0, i])
sigma_z_y2 = 0
for i in range(n): sigma_z_y2 += float(Z[0, i] * Y[0, i] * Y[0, i])
sigma_z_x_y = 0
for i in range(n): sigma_z_x_y += float(Z[0, i] * X[0, i] * Y[0, i])
sigma_z_x = 0
for i in range(n): sigma_z_x += float(Z[0, i] * X[0, i])
sigma_z_y = 0
for i in range(n): sigma_z_y += float(Z[0, i] * Y[0, 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)
return res
def matching_3D(X, Y, Z):
n = np.size(X)
print(n)
res = get_res(X, Y, Z, n)
# 输出方程形式（注意精度，否则系数无法完全显示！！！！！）
print("z=%.32s*x^2%.32s*xy%.32s*y^2%.32s*x%.32s*y%.32s" % (
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')
plt.show()
def test():
#验证准确性
X, Y, Z = read_csv()
x=X[15,5]
y=Y[15,5]
print(x)
print(y)
z = -4.437e-14 * pow(x,2) + 6.717e-10 * x*y + 2.5213e-11 * pow(y,2) + 1.477e-8 * x - 5.943e-8 * y + 2.129e-5
print(z)
if __name__ == '__main__':
X,Y,Z = read_csv()
X = np.reshape(X, (1, -1))
Y = np.reshape(Y, (1, -1))
Z = np.reshape(Z, (1, -1))
# print("X",X)
# print("Y", Y)
# print("Z", Z)
matching_3D(X, Y, Z)
test()
python结果
z=-4.437671280867266e-14*x^2+6.717093852298463e-10*xy+2.5212810541556962e-11*y^2+1.4772941549929638e-08*x-5.943476386001641e-08*y+2.1293990002750532e-05
参考文章：
https://blog.csdn.net/Haipai1998/article/details/85345823
https://blog.csdn.net/weixin_44112790/article/details/89287521
扫一扫
2597
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
