全面地介绍将二维数据矩阵插值拟合为二维曲面方程—

这篇具有很好参考价值的文章主要介绍了全面地介绍将二维数据矩阵插值拟合为二维曲面方程——基于matlab。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

前言

要利用matlab将二维数据矩阵拟合成一个二维曲面方程，首先，我们得明白什么是插值和拟合。下面，从对洋流数据的处理入手，先以一个基本的流程为例，进行问题的概述，然后介绍两种具体的方法。

matlab根据数据拟合曲面,研究生学习,matlab,数学建模,学习,考研,经验分享,学习方法,矩阵

原理初析

插值

根据已知的数据点，生成一个连续的函数，可以在任意位置求值。matlab 提供了多种插值方法，如线性插值、 $B$ -样条插值、多项式插值等。

这些方法的优缺点和适用场合各有不同，一般来说，线性插值比较简单和快速，但可能不够平滑； $B$ -样条插值比较平滑和灵活，但可能不够精确；多项式插值比较精确和光滑，但可能出现奇异或震荡现象。具体选择哪种方法，要根据数据的特点和目的来决定。

matlab中，可以使用 $in t er p 2$ 函数或 $g r i dd e d I n t er p o l an t$ 类来进行二维插值，它们都可以指定不同的插值方法。例如，如果已知的数据矩阵是 ${V}$ ，想要在 ${Xq}$ 和 ${Yq}$ 的网格上进行插值，可以使用以下代码：


% 线性插值
Vq1 = interp2(V,Xq,Yq,'linear');
% B-样条插值
Vq2 = interp2(V,Xq,Yq,'spline');
% 多项式插值
Vq3 = interp2(V,Xq,Yq,'pchip');

拟合

根据插值得到的函数，寻找一个合适的数学模型，可以用一个简洁的表达式来近似描述数据的规律。

matlab提供了多种拟合方法，如线性回归、非线性回归、最小二乘法等。这些方法的原理和应用也各有不同，一般来说，线性回归比较简单和稳定，但可能不够灵活；非线性回归比较灵活和通用，但可能不够稳定；最小二乘法比较优化和精确，但可能不够简洁。具体选择哪种方法，要根据数据的分布和形状来决定。matlab中，可以使用 $f i t$ 函数或 $f i tt y p e$ 类来进行二维拟合，它们都可以指定不同的拟合方法。例如，如果插值得到的函数是 ${Vq}$ ，想要用一个二次多项式来拟合，可以使用以下代码：

% 线性回归
ft1 = fittype('poly22'); % 二次多项式
f1 = fit([Xq(:),Yq(:)],Vq1(:),ft1); % 拟合
% 非线性回归
ft2 = fittype('a*exp(b*x)*sin(c*y)'); % 非线性函数
f2 = fit([Xq(:),Yq(:)],Vq2(:),ft2); % 拟合
% 最小二乘法
ft3 = fittype('a*x^2+b*y^2+c*x*y+d*x+e*y+f','problem','f'); % 二次多项式，f是固定参数
f3 = fit([Xq(:),Yq(:)],Vq3(:),ft3,'problem',0); % 拟合

拟合完成后，可以用 f1，f2，f3等变量来表示拟合得到的二维曲面方程，也可以用 $coe ff v a l u es$ 函数来提取方程的系数。

例如，如果想要得到 f1的表达式，可以使用以下代码：

% 提取系数
c1 = coeffvalues(f1);
% 构造表达式
s1 = sprintf('V = %.4f*X^2 + %.4f*Y^2 + %.4f*X*Y + %.4f*X + %.4f*Y + %.4f',c1(1),c1(2),c1(3),c1(4),c1(5),c1(6));

拟合优度的评估

对于方程的拟合，如果要做可行性分析的话，那必然要对拟合优度进行评估。

拟合优度是一种衡量拟合结果与实际数据之间的差异的指标。拟合优度越高，说明拟合结果越能反映数据的真实情况。拟合优度的评估方法有很多种，常见的有残差分析、决定系数、均方根误差等。不同的评估方法有不同的含义和应用场景，比如残差分析可以检验拟合函数的合理性，决定系数可以衡量拟合函数对数据变化的解释程度，均方根误差可以衡量拟合函数的预测精度等。

具体实例的展开

下面，将以两种方法进行具体问题的求解。

B-样条插值+最小二乘法

B-样条插值法

$B$ -样条（B-spline）是一种基于多项式的插值方法，它可以构造出平滑且连续的曲线或曲面。b-样条的优点是可以灵活地控制插值的阶数、节点、边界条件等，从而得到不同的插值效果。b-样条的缺点是计算量较大，且难以表示复杂的形状。

$B$ -样条的数学定义如下：

设有 $n + 1$ 个数据点 $\left(x_i,y_i\right)$ ，其中 $i=0,1,\dots,n$ ，且 $x_0<x_1<\dots<x_n$ 。 $B$ -样条插值法的目标是找到一个函数 $f\left(x\right)$ ，使得 $f\left(x_i\right)=y_i$ ，对所有 $i$ 成立，且 $f\left(x\right)$ 在每个区间 $\left[x_i,x_{i+1}\right]$ 上是一个 $k$ 阶多项式，其中 $k$ 是给定的插值阶数，通常为 $3$ 或 $4$ 。

为了构造这样的函数，我们需要引入一个节点序列 $t_0,t_1,\dots,t_{n+k}$ ，其中 $t_i$ 是一个实数，且满足 $t_0\le t_1\le \dots\le t_{n+k}$ 。节点序列的作用是划分数据点的区间，并确定每个区间上的多项式的系数。根据节点序列，我们可以定义一组基函数 $B_{i,k}(x)$ ，其中 $i=0,1,\dots,n，k=1,2,\dots$ ，满足以下递推公式：

$B_{i,1}(x)=\begin{cases}1,&t_i\le x<t_{i+1}\\0,&\text{otherwise}\end{cases}$
$B_{i,k}(x)=\frac{x-t_i}{t_{i+k-1}-t_i}B_{i,k-1}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}B_{i+1,k-1}(x)$

基函数 $B_{i,k}(x)$ 的意义是在节点序列 $t_i,t_{i+1},\dots,$ 上的 $k$ 阶 $B$ -样条，它是一个非负的、分段的、最高为 $k - 1$ 次的多项式，且在 $t_i,t_{i+1}]$ 上为 $1$ ，在 $t_{i+k-1},t_{i+k}]$ 上为 $0$ ，在其他区间上线性变化。

基于基函数，我们可以将插值函数 $f (x)$ 表示为：

$f(x)=\sum_{i=0}^n a_i B_{i,k}(x)$
其中 $a_i$ 是待求的系数，它们可以通过解线性方程组 $\sum_{i=0}^n a_i B_{i,k}(x_j)=y_j$

对 $j=0,1,\dots,n$ 得到。

这样，我们就得到了一个 $B$ -样条插值函数 $f (x)$ ，它可以在任意给定的 $x$ 值处计算出对应的 $y$ 值。

如果我们要得到一个二维曲面方程 $f (x, y)$ ，我们可以先沿着 $x$ 方向对每个 $y$ 值进行 $B$ -样条插值，得到一组函数

$f_i(y)=\sum_{j=0}^n a_{ij} B_{j,k}(x)$

其中 $i=0,1,\dots,m，$ 然后再沿着 $y$ 方向对每个 $x$ 值进行 $B$ -样条插值，得到最终的函数

$f(x,y)=\sum_{i=0}^m f_i(y) B_{i,k}(y)$

这样，我们就得到了一个双向的 $B$ -样条插值曲面，它可以在任意给定的 $(x, y)$ 值处计算出对应的 $z$ 值。

lsqcurvefit 拟合法

matlab根据数据拟合曲面,研究生学习,matlab,数学建模,学习,考研,经验分享,学习方法,矩阵

接下来，介绍一种主要的最小二乘法的拟合算法在matlab里的应用—— $l s q c u r v e f i t$ 拟合法。

$l s q c u r v e f i t$ 是 MATLAB 中的一个函数，它可以用来对数据进行非线性曲线拟合。 $l s q c u r v e f i t$ 的优点是可以灵活地指定拟合函数的形式，且可以考虑数据的误差和噪声。 $l s q c u r v e f i t$ 的缺点是需要提供拟合函数的初始参数，且可能存在局部最优解的问题。

设有 $n$ 个数据点 $x_i,y_i)$ ，其中 $i=1,2,\dots,n$ ，且 $x_i$ 和 $y_i$ 可以是标量或向量。 $l s q c u r v e f i t$ 的目标是找到一个函数 $f (x, a)$ ，使得

$\sum_{i=1}^n \|y_i-f(x_i,a)\|^2$

最小，其中 $a$ 是待求的参数向量，它可以有一些约束条件，比如上下界、线性或非线性不等式或等式。

为了求解这样的问题， $l s q c u r v e f i t$ 使用了一种基于信赖域的迭代算法。

它的基本思想是在每一步，以当前的参数值为中心，构造一个信赖域，然后在信赖域内用一个二次模型来近似目标函数，求解该二次模型的最优解，作为下一步的参数值。这个过程重复进行，直到满足收敛条件或达到最大迭代次数。信赖域的大小会根据每一步的优化效果进行调整，以保证算法的稳定性和效率。

$l s q c u r v e f i t$ 在具体编程的时候需要提供以下输入：

fun：一个函数句柄，表示拟合函数 $f (x, a)$ ，它接受两个输入参数，分别是自变量 $x$ 和参数向量 $a$ ，并返回一个输出值，表示因变量 $y$ 。
x0：一个向量，表示参数向量 $a$ 的初始值，它决定了算法的起始点，对于非线性问题，不同的初始值可能导致不同的最优解，所以应该尽量选择一个合理的初始值。
xdata：一个向量或矩阵，表示自变量 $x$ 的数据，它的大小应该与因变量 $y$ 的数据相匹配。
ydata：一个向量或矩阵，表示因变量 $y$ 的数据，它的大小应该与自变量 $x$ 的数据相匹配。
lb：一个向量，表示参数向量 $a$ 的下界，它的大小应该与参数向量 $a$ 相同，如果没有下界，可以用 Inf 表示。
ub：一个向量，表示参数向量 $a$ 的上界，它的大小应该与参数向量 $a$ 相同，如果没有上界，可以用 Inf 表示。

$l s q c u r v e f i t$ 的输出是一个向量，表示参数向量 $a$ 的最优值，它使得目标函数的平方和达到最小。

如果我们要得到一个二维曲面方程 $f (x, y, a)$ ，我们可以先将自变量 $x$ 和 $y$ 组合成一个矩阵，然后将其作为 xdata 输入 $l s q c u r v e f i t$ ，得到参数向量 $a$ 的最优值，然后将其代入拟合函数 $f (x, y, a)$ ，得到最终的函数 $f (x, y)$ 。这样，我们就得到了一个非线性的曲面拟合方程，它可以在任意给定的 $(x, y)$ 值处计算出对应的 $z$ 值。

数学实现

为了演示插值和拟合的数学实现，使用一个示例数据矩阵，它包含了 $x$ ， $y$ 位置方向的， $u$ ， $v$ 方向速度的洋流数据，大小为 $50\times 50$ 。

这里，假设你的洋流数据是从官方网站下载的，所以我的这个算法就假设你直接 load 这个洋流数据文件。

读取数据矩阵，并提取 $x ， y ， u ， v$ 四个变量，分别表示位置和速度。这里我假设数据矩阵已经保存在一个名为 data.mat 的文件中，各位小伙伴们可以根据自己的实际情况修改文件名和路径哒。

load data.mat
x = data(:,1);
y = data(:,2);
u = data(:,3);
v = data(:,4);

将 $x ， y ， u ， v$ 四个变量转换为矩阵形式，以便进行插值和拟合。这里我假设 $x ， y$ 是均匀分布的，当然也可以根据自己的实际情况修改矩阵的大小和范围。

X = reshape(x,50,50);
Y = reshape(y,50,50);
U = reshape(u,50,50);
V = reshape(v,50,50);

计算 $z$ 变量，表示速度的大小，它是 $u ， v$ 的平方和的开方。这个变量将作为插值和拟合的因变量，表示洋流的强度。

Z = sqrt(U.^2+V.^2);

使用 $B$ -样条插值法，对 $z$ 变量进行插值，得到一个二维曲面方程 $f_1(x,y)$ 。

这里，我使用 MATLAB 的 $C u r v e F i tt in g T oo l b o x$ 中的 $f i t$ 函数，指定插值方法为 $‘ biha r m o ni c in t er p ’$ ，表示双向的三次 $B$ -样条插值。也可以使用其他的插值方法，比如 $‘ t hin pl a t e in t er p ’$ ，表示双向的二次 $B$ -样条插值，或者 $‘ l in e a r in t er p ’$ ，表示双向的线性插值。

$f i t$ 函数的输出是一个 $c f i t$ 对象，它包含了插值函数的信息，也可以使用 $f e v a l$ 函数来计算插值函数在任意 $(x, y)$ 值处的 $z$ 值，或者使用 $f or m u l a$ 函数来查看插值函数的具体表达式。

注意，由于 b-样条插值函数是分段的，它的表达式可能很长，且难以理解，所以我不建议各位友友直接查看它的表达式，而是使用 $f e v a l$ 函数来计算它的值。

f1 = fit([X(:),Y(:)],Z(:),'biharmonicinterp');
Z1 = feval(f1,X,Y);

接下来，使用 $l s q c u r v e f i t$ 拟合法，对 $z$ 变量进行拟合，得到一个二维曲面方程 $f_2(x,y,a)$ 。这里我假设拟合函数的形式为 $f_2(x,y,a)=a_1\sin(a_2x+a_3y)+a_4\cos(a_5x+a_6y)+a_7$

，其中 $a=(a_1,a_2,a_3,a_4,a_5,a_6,a_7)$ 是待求的参数向量。

当然，你也可以根据自己的实际情况修改拟合函数的形式，但请注意，拟合函数应该是非线性的，且能够较好地拟合数据的特征。

$l s q c u r v e f i t$ 函数的输入包括拟合函数的函数句柄，参数向量的初始值，自变量和因变量的数据，参数向量的上下界等。 $l s q c u r v e f i t$ 函数的输出是参数向量的最优值，它使得目标函数的平方和达到最小。

根据参数向量的最优值，我们可以得到拟合函数的具体表达式，也可以使用 $f e v a l$ 函数来计算拟合函数在任意 $(x, y)$ 值处的 $z$ 值。

fun = @(a,xdata) a(1)*sin(a(2)*xdata(:,1)+a(3)*xdata(:,2))+a(4)*cos(a(5)*xdata(:,1)+a(6)*xdata(:,2))+a(7);
x0 = [1,1,1,1,1,1,1];
xdata = [X(:),Y(:)];
ydata = Z(:);
lb = [-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf];
ub = [Inf,Inf,Inf,Inf,Inf,Inf,Inf];
a = lsqcurvefit(fun,x0,xdata,ydata,lb,ub);
Z2 = feval(fun,a,xdata);
Z2 = reshape(Z2,50,50);
f2 = @(x,y) a(1)*sin(a(2)*x+a(3)*y)+a(4)*cos(a(5)*x+a(6)*y)+a(7);

代码

接下来，是一个整个的代码实现。

% Load data matrix and extract variables
load data.mat
x = data(:,1);
y = data(:,2);
u = data(:,3);
v = data(:,4);

% Convert variables to matrix form
X = reshape(x,50,50);
Y = reshape(y,50,50);
U = reshape(u,50,50);
V = reshape(v,50,50);



% Calculate z variable as the magnitude of velocity 
Z = sqrt(U^2+V^2);
% Use b-spline interpolation method to interpolate z variable and get a surface equation f1(x,y) 
f1 = fit([X(:),Y(:)],Z(:),'biharmonicinterp'); Z1 = feval(f1,X,Y);
% Use lsqcurvefit fitting method to fit z variable and get a surface equation f2(x,y,a) 
fun = @(a,xdata) a(1)*sin(a(2)*xdata(:,1)+a(3)*xdata(:,2))+a(4)*cos(a(5)*xdata(:,1)+a(6)*xdata(:,2))+a(7); 
x0 = [1,1,1,1,1,1,1]; 
xdata = [X(:),Y(:)]; 
ydata = Z(:); 
lb = [-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf]; 
ub = [Inf,Inf,Inf,Inf,Inf,Inf,Inf]; 
a = lsqcurvefit(fun,x0,xdata,ydata,lb,ub); 
Z2 = feval(fun,a,xdata); 
Z2 = reshape(Z2,50,50); 
f2 = @(x,y) a(1)*sin(a(2)*x+a(3)*y)+a(4)*cos(a(5)*x+a(6)*y)+a(7);
% Plot the original data, the interpolated surface, and the fitted surface 
figure 
subplot(1,3,1) 
surf(X,Y,Z) 
title('Original Data') 
xlabel('x'),ylabel('y'),zlabel('z') 
subplot(1,3,2) 
surf(X,Y,Z1) 
title('Interpolated Surface') 
xlabel('x'), ylabel('y'),zlabel('z') 
subplot(1,3,3) 
surf(X,Y,Z2) 
title('Fitted Surface') 
xlabel('x'),ylabel('y'),zlabel('z')
% Compare the interpolation and fitting results with the original data % Calculate the mean squared error (MSE) and the coefficient of determination (R-squared) for each method 
mse1 = mean((Z(:)-Z1(:)).^2); 
mse2 = mean((Z(:)-Z2(:)).^2); 
r1 = 1 - sum((Z(:)-Z1(:)).^2)/sum((Z(:)-mean(Z(:))).^2); 
r2 = 1 - sum((Z(:)-Z2(:))^.2)/sum((Z(:)-mean(Z(:)))^.2); 
fprintf('The MSE and R-squared for b-spline interpolation are: %.4f and %.4f\n',mse1,r1); 
fprintf('The MSE and R-squared for lsqcurvefit fitting are: %.4f and %.4f\n',mse2,r2);
% Display the interpolation and fitting functions 
fprintf('The interpolation function is: f1(x,y) = %s\n',formula(f1));
fprintf('The fitting function is: f2(x,y) = %s\n',formula(f2));

B-样条插值+非线性回归法

这个方法就不着重讲原理分析了，讲一下拟合评估怎么做，还附带讲一讲洋流的可视化。

数据的预处理

同样，假设你已经有了一个 $50 \times 50$ 的二维数据矩阵，比如洋流数据矩阵，其中每个元素表示某个位置的洋流速度。为了进行插值和拟合，你需要将这个矩阵转换为三个向量，分别表示 $x$ 坐标、 $y$ 坐标和 $z$ 值。可以使用 matlab 的 $res ha p e$ 函数来实现这一转换:

% 假设 Z 是 50 x 50 的洋流数据矩阵
n = 50; % 矩阵的大小
x = repmat(1:n,1,n); % 重复 1 到 n，共 n 次，得到 x 坐标向量
y = reshape(repmat(1:n,n,1),1,n^2); % 重复 1 到 n，每个 n 次，然后转换为行向量，得到 y 坐标向量
z = reshape(Z,1,n^2); % 将 Z 转换为行向量，得到 z 值向量

% 假设 x, y, z 是已经转换好的三个向量，表示数据点的坐标和值
m = 100; % 网格的大小
xi = linspace(min(x),max(x),m); % 生成 x 方向的网格坐标
yi = linspace(min(y),max(y),m); % 生成 y 方向的网格坐标
[XI,YI] = meshgrid(xi,yi); % 生成网格矩阵
ZI = griddata(x,y,z,XI,YI,'cubic'); % 使用样条插值得到网格上的值

非线性回归的实现

为了使用 $f i t$ 函数进行非线性回归，你需要指定拟合函数的形式和参数的初始值或范围。拟合函数的形式可以是一个匿名函数，也可以是一个字符串，表示函数的表达式。参数的初始值或范围可以是一个向量，也可以是一个矩阵，表示参数的下限和上限：

% 假设 x, y, z 是已经转换好的三个向量，表示数据点的坐标和值
f = @(a,b,c,d,x,y) a*x.^2 + b*y.^2 + c*x.*y + d; % 定义拟合函数的形式，这里是一个二次曲面
s = [1 1 1 1]; % 定义参数的初始值
f = fit([x,y],z,f,'StartPoint',s); % 使用非线性回归得到拟合结果

拟合优度的评估

前面已提及了拟合优度评估的原理，这里讲讲代码的实现。为了使用 $f i t$ 函数的输出来评估拟合优度，你可以查看拟合结果的属性，其中包含了一些拟合优度的指标：

% 假设 f 是拟合结果
R2 = f.rsquare; % 得到决定系数
RMSE = f.rmse; % 得到均方根误差
disp(['R-squared = ',num2str(R2)]); % 显示决定系数
disp(['RMSE = ',num2str(RMSE)]); % 显示均方根误差

曲面方程的求解

为了使用 $f i t$ 函数的输出来求解曲面方程，你可以使用 $f$ 函数来计算任意位置的 $z$ 值，或者使用 $f e v a l$ 函数来计算网格上的 $z$ 值:

% 假设 f 是拟合结果，XI, YI 是网格矩阵
z1 = f(10,20); % 计算 x = 10, y = 20 时的 z 值
z2 = feval(f,XI,YI); % 计算网格上的 z 值
disp(['z(10,20) = ',num2str(z1)]); % 显示 z(10,20) 的值

矢量场的绘制

对于洋流的可视化，因为洋流由 $x, y, u, v$ 这四组向量组成，其中前两个定位洋流的每个坐标点的位置，后两个代表每个坐标点水平方向和竖直方向的速度，所以可以用矢量场来可视化，为了绘制矢量场，可以使用 $q u i v er$ 函数来绘制箭头，表示矢量的方向和大小。

在绘制的时候，需要指定箭头的起点和终点的坐标，以及箭头的颜色、线宽、比例等属性：

% 假设 XI, YI, ZI 是网格矩阵，表示插值后的数据
[U,V,W] = gradient(ZI); % 计算网格上的梯度，表示矢量的方向和大小
quiver3(XI,YI,ZI,U,V,W,'r','LineWidth',2,'AutoScaleFactor',0.5); % 绘制矢量场，使用红色的箭头，线宽为 2 ，比例因子为 0.5
xlabel('x'); % 标注 x 轴
ylabel('y'); % 标注 y 轴
zlabel('z'); % 标注 z 轴
title('Vector field of the fitted surface'); % 标注标题

尾声

本文由一个通用方法进行问题的介绍，然后，使用了两种组合方法进行具体问题求解的展开。

将二维数据矩阵转化为二维曲面方程有很多方法，可能本文的解答有局限，请多多见谅。

最后，祝君功不唐捐，玉汝于成！

matlab根据数据拟合曲面,研究生学习,matlab,数学建模,学习,考研,经验分享,学习方法,矩阵文章来源地址https://www.toymoban.com/news/detail-847975.html

参考

MathWorks. (2021). Interpolate scattered data.
MathWorks. (2021). Representing a matrix as a surface.
Wikipedia contributors. (2024, March 2). Bilinear interpolation. In Wikipedia, The Free Encyclopedia. Retrieved 02:41, March 4, 2024.
Kiusalaas, J. (2013). Interpolation in Two and Three Dimensions. In Numerical Methods in Engineering with Python 3 (pp. 97-122). Cambridge: Cambridge University Press. doi:10.1017/CBO9781139108278.006
de Boor, C. (1978). Biharmonic spline interpolation. In A Practical Guide to Splines (pp. 211-247). New York: Springer-Verlag.
MathWorks. (n.d.). Two-Dimensional Interpolation. Retrieved March 4, 2024.
XlXtrFun. (n.d.). Excel 2D interpolation function. Retrieved March 4, 2024.

到了这里，关于全面地介绍将二维数据矩阵插值拟合为二维曲面方程——基于matlab的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！