clc;clear;close all; pt = pcread('11-左水平尾翼上表面.pcd'); P = pt.Location; x = P(:,1); y = P(:,2); z = P(:,3); % 定义LSCM矩阵 L = [sum(x.^2),sum(x.*y),sum(x); sum(x.*y),sum(y.^2),sum(y); sum(x),sum(y),length(x)]; % 定义右侧向量 b = [sum(x.*z);sum(y.*z);sum(z)]; % 计算最小二乘解 a = inv(L)*b; % 定义曲面函数 f = @(x,y) a(1)*x.^2 + a(2)*x.*y + a(3)*x + a(4)*y.^2 + a(5)*y + a(6); % 绘制数据点和曲面 xx = min(x) :1 : max(x); yy = max(y) :-1: min(y); [X,Y] = meshgrid(xx, yy) ; Z = f(X,Y); scatter3(x,y,z,'filled') hold on surf(X,Y,Z) xlabel('X') ylabel('Y') zlabel('Z')。怎么改正
时间: 2023-10-22 13:03:20 浏览: 78
MATLAB经验总结-2014-7-18(持续更新).pdf
这段代码的问题是在最小二乘解的计算中使用了求逆运算,这样做不仅计算量大,而且容易出现数值不稳定的情况。可以使用奇异值分解(SVD)的方法来求解最小二乘解,代码如下:
clc;clear;close all;
pt = pcread('11-左水平尾翼上表面.pcd');
P = pt.Location;
x = P(:,1);
y = P(:,2);
z = P(:,3);
% 定义LSCM矩阵
L = [x.^2, x.*y, x, y.^2, y, ones(size(x))];
% 定义右侧向量
b = z;
% 计算最小二乘解
[U,S,V] = svd(L,'econ');
a = V * diag(1./diag(S)) * U' * b';
% 定义曲面函数
f = @(x,y) a(1)*x.^2 + a(2)*x.*y + a(3)*x + a(4)*y.^2 + a(5)*y + a(6);
% 绘制数据点和曲面
xx = min(x) :1 : max(x);
yy = max(y) :-1: min(y);
[X,Y] = meshgrid(xx, yy) ;
Z = f(X,Y);
scatter3(x,y,z,'filled')
hold on
surf(X,Y,Z)
xlabel('X')
ylabel('Y')
zlabel('Z')
这样可以避免求逆操作,计算更稳定,也更高效。
阅读全文