用python拟合三维数据的曲面，并画出曲面的等高线图，求出曲面的波峰和波谷

主要用 scipy.optimize.curve_fit 拟合曲面
用matplotlib画图

试试这段逻辑，先定义一个二次曲面，然后用curve_fit拟合，其他里面有备注，试试看


import numpy as np
from scipy.optimize import curve_fit
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# 定义二次曲面模型
def model(data, a, b, c, d, e, f):
    x, y = data
    return a * x**2 + b * y**2 + c * x * y + d * x + e * y + f

# 生成一些示例数据
x = np.random.rand(100) * 10
y = np.random.rand(100) * 10
z = model((x, y), 0.5, 0.3, 0.1, 0.2, 0.1, 1.0) + np.random.rand(100) - 0.5

# 拟合曲面
params, _ = curve_fit(model, (x, y), z)

# 生成曲面和等高线图的数据
x_range = np.linspace(x.min(), x.max(), 100)
y_range = np.linspace(y.min(), y.max(), 100)
x_grid, y_grid = np.meshgrid(x_range, y_range)
z_grid = model((x_grid, y_grid), *params)

# 绘制曲面图
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='b')  # 原始数据
ax.plot_surface(x_grid, y_grid, z_grid, color='r', alpha=0.3)  # 拟合曲面
plt.show()

# 绘制等高线图
plt.contour(x_grid, y_grid, z_grid, levels=20)
plt.show()

# 通过求导找出极值点
# 请注意，在这里我们假设模型是二次曲面，如果使用其他模型，求导的过程将会不同
# 对于二次曲面，极值点就是曲面的顶点，我们可以直接通过公式计算出来
x_peak = -(params[3] * params[1] - params[2] * params[4]) / (4 * params[0] * params[1] - params[2]**2)
y_peak = -(params[4] * params[0] - params[2] * params[3]) / (4 * params[0] * params[1] - params[2]**2)
z_peak = model((x_peak, y_peak), *params)
print(f'Peak: x={x_peak}, y={y_peak}, z={z_peak}')

如果有帮助，请点击一下采纳该答案～谢谢

你可以看下这个问题的回答https://ask.csdn.net/questions/7799516
这篇博客也不错, 你可以看下Python绘制曲面和曲面对应的等高线图
除此之外, 这篇博客: 利用Python提取函数图像数据并拟合曲线中的 3. 曲线拟合 部分也许能够解决你的问题, 你可以仔细阅读以下内容或跳转源博客中阅读:

import numpy as np
import pandas as pd
from sympy import *
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%config ZMQInteractiveShell.ast_node_interactivity = "all"

df = pd.read_excel('D:/learning/Function coordinates/test/p1_data.xlsx', index_col = 0)
X = df['x'].values
Y = df['y'].values

def function1(x,a1,a2,a3,b1,b2,b3,c11,c12,c2,c3):
    return 2*a1 / (np.exp(-(x-b1)/c11)+np.exp((x-b1)/c12)) +\
            a2 * (b2**2) / (x**2) * np.exp(-((x-b2)/(x*c2))**2) +\
            a3 * (b3**2) / (x**2) * np.exp(-((x-b3)/(x*c3))**2)
popt, pcov = curve_fit(function1, X, Y, bounds=([0.,0.,0.,380.,400.,400.,0.,0.,0.,0.], 
                                            [1.,0.6,0.6,520.,800.,800.,20.,100.,100.,100.]))
print(popt)
# print(pcov)

Y1 = 2*popt[0] / (np.exp(-(X-popt[3])/popt[6])+np.exp((X-popt[3])/popt[7]))
Y2 = popt[1] * (popt[4]**2) / (X**2) * np.exp(-((X-popt[4])/(X*popt[8]))**2)
Y3 = popt[2] * (popt[5]**2) / (X**2) * np.exp(-((X-popt[5])/(X*popt[9]))**2)

plt.plot(X, Y, label='data')
plt.plot(X, function1(X,*popt), label='fit')
plt.plot(X, Y1, label='T1')
plt.plot(X, Y2, label='T2')
plt.plot(X, Y3, label='T3')
plt.legend(loc='upper right')

p11
拟合曲线里面具体的计算原理不太了解。
拟合时 curve_fit() 中 bounds 里面的两个列表对应系数 a1,a2,a3,b1,b2,b3,c11,c12,c2,c3 的最小值和最大值，其数值会在一定程度影响拟合效果，打印的 popt 是系数的值。
拟合公式为：
在这里插入图片描述
顺便计算了一下拟合曲线的绝对误差：

In: abs(Y - function1(X,*popt)).sum()
Out: 16.800576957203255

还做了个用6个高斯拟合的：
在这里插入图片描述

def gs6(x,a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,c1,c2,c3,c4,c5,c6):
    return a1*np.exp(-((x-b1)/c1)**2)+\
            a2*np.exp(-((x-b2)/c2)**2)+\
            a3*np.exp(-((x-b3)/c3)**2)+\
            a4*np.exp(-((x-b4)/c4)**2)+\
            a5*np.exp(-((x-b5)/c5)**2)+\
            a6*np.exp(-((x-b6)/c6)**2)
popt, pcov = curve_fit(gs6, X, Y, bounds=([0.,0.,0.,0.,0.,0.,400.,400.,400.,400.,400.,400.,0.,0.,0.,0.,0.,0.], 
                                            [1.2,1.2,1.2,0.6,0.6,0.6,510.,510.,510.,800.,800.,800.,20.,20.,20.,100.,90.,90.]))
print(popt)
# print(pcov)

Y1 = popt[0]*np.exp(-((X-popt[6])/popt[12])**2)
Y2 = popt[1]*np.exp(-((X-popt[7])/popt[13])**2)
Y3 = popt[2]*np.exp(-((X-popt[8])/popt[14])**2)
Y4 = popt[3]*np.exp(-((X-popt[9])/popt[15])**2)
Y5 = popt[4]*np.exp(-((X-popt[10])/popt[16])**2)
Y6 = popt[5]*np.exp(-((X-popt[11])/popt[17])**2)

plt.plot(X, Y, label='data')
plt.plot(X, gs6(X,*popt), label='fit')
plt.plot(X, Y1, label='T1')
plt.plot(X, Y2, label='T2')
plt.plot(X, Y3, label='T3')
plt.plot(X, Y4, label='T4')
plt.plot(X, Y5, label='T5')
plt.plot(X, Y6, label='T6')
plt.legend(loc='upper right')

在这里插入图片描述

In: abs(Y - gs6(X,*popt)).sum()
Out: 6.379130579724297