panlei 发表于 2009-6-25 19:12

求助-曲面拟合

function ff=main(x,p,y,q,z,xx,yy)
% x y z 坐标向量 长度要一样
% p ,q 为拟合函数中x,y 的系数
% xx yy 为 需要拟合的数据 给出(x,y) 坐标 求z
x=input('x=');
p=input('p=');
y=input('y=');
q=input('q=');
z=input('z=');
xx=input('xx=');
yy=input('yy=');
A=leftmatrix (x,p,y,q);             % A*a_n=B
B=rightmatrix(x,p,y,q,z);
%a_n=inv(A)*B;
a_n=A\B;                        % 求a_n    inv(A)*B 效果不好 (存疑)
for i=1 : length(a_n)             % 把长为p*q 的 a_n 列向量 转还 成p x q的 aa 矩阵
    ii=quotient(i-1,q)+1;         % quotient求商
    jj=mod(i-1,q)+1; aa(ii,jj)=a_n(i,1);
end
                                 
ff=0;                              % ff 是 xx,yy 带入所拟合的函数 求出 z
for i=1 : p                        % 函数为 aa(i,j)*x^i*y^j(i=0...p,j=0...q)
    for j=1 : q                  % aa 为系数p x q 的矩阵
      ff=ff+aa(i,j) * xx.^(i-1) * yy.^(j-1);
    end
end
function A=leftmatrix(x,p,y,q)
% A*a=Ba 即为系数列矩阵
% A为左边(p-1)(q-1) 乘 (p-1)(q-1) 的矩阵
% x,y 为长度一样的列矩阵 也就是给定离散点的x,y坐标
% p,q为拟合的函数中x,y的指数
%x=input('x=');
%y=input('y=');
m=length(x);
if ((nargin~=4) & (m~=length(y))), error('error check check!'); end
A_length=p*q;                        % A 为p*q阶的方阵
A=zeros(A_length,A_length);          % 赋值0
for i=1 : p*q
    for j= 1 : p*q
      x_z=quotient(j-1,q)+quotient(i-1,q);   % x 的指数   quotient为求商
      y_z=mod(j-1,q)+mod(i-1,q);            % y 的指数
      A(i,j)=qiuhe(x,x_z,y,y_z);            
    end
end
function he=qiuhe(x,p,y,q,z)
% he    x^p*y^q 从1->m的和
% x,y 行向量 长度相同; p,q 为x,y系数
%x=; y=[ 3 4]; p=2; q= 2;
m=length(x);
if (nargin<4 )&(m~=length(y))            %输入量至少为四,x,y行向量长度必需一样
    error('error check check!');
end
if nargin==4, z=ones(m,1); end%没有 z , 默认为单位行向量
he=0;
for i=1:m
    he=he+x(i).^p * y(i).^q*z(i);   % 1-->m 求和
end
function sh=quotient(x,y)
% sh 为 x/y 的商
sh=(x-mod(x,y))/y;
function B=rightmatrix(x,p,y,q,z)
% A*a=B
% B为一个列向量 长为p*q
% x y z 为点的坐标 , p q 为x y指数
if nargin~=5,error('error check check! rightmatrix'); end
B=zeros(p*q,1);
for i=1 : p*q
    x_z=quotient(i-1,q); y_z=mod(i-1,q); B(i,1)=qiuhe(x,x_z,y,y_z,z);
end

以上是程序,哪位大侠帮帮我看看哪里有问题,输入数据的时候总是出现问题
??? Error using ==> mtimes
Inner matrix dimensions must agree.
Error in ==> qumiannihe at 27
      ff=ff+aa(i,j) * xx.^(i-1) * yy.^(j-1);
在此先谢过了!我的邮箱panlei00@126.com

ChaChing 发表于 2009-6-25 19:56

Ref : 5F
常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html

还有这帖怎好像旧帖!?
http://forum.vibunion.com/forum/thread-43640-1-1.html

[ 本帖最后由 ChaChing 于 2009-6-25 20:42 编辑 ]

panlei 发表于 2009-6-25 20:50

感谢你的信息,不过我们出现的问题不同呀,本人这方面较差,还望仔细说明,不胜感激!

panlei 发表于 2009-6-26 11:07

感谢你的信息,不过我们出现的问题不同呀,本人这方面较差,还望仔细说明,不胜感激!
我的输入
x=; p=4; y=; q=4;
z=; xx=; yy=;
Warning: Matrix is singular to working precision.
> In qumiannihe at 17
??? Error using ==> mpower
Matrix must be square.
Error in ==> qumiannihe at 27
      ff=ff+aa(i,j) * xx^(i-1) * yy^(j-1);
请高手帮帮忙,先谢过了!

[ 本帖最后由 ChaChing 于 2009-6-26 22:26 编辑 ]

ChaChing 发表于 2009-6-26 22:37

说真的, 这种很长一大串的程序, 个人很懒得看! 没法, 时间有限!
1.qumiannihe是什麽? 没交代! 当我是半仙, 用算的?:@)
2."qumiannihe at 27", 主程序那里来的27行?
3.怎两次的报错会不同?
...

[ 本帖最后由 ChaChing 于 2009-6-26 22:51 编辑 ]

panlei 发表于 2009-6-27 12:26

原理:
给定一组坐标(x,y,z) 拟合成函数
(x,y,z 为行向量)
   (aa 为p x q 的系数矩阵)
由最小二乘法得:
点( ….. )是多元函数


S( ….. )=
(w(x,y)为权函数,默认为1)
的极小点,从而 ….. 必须满足方程组
                  (ij=0…pq)
对S函数求偏导,移项之后得

其中 , ,
也就是A*a_n=B

a_n=

B=

A 为pq 阶矩阵

function A=leftmatrix(x,p,y,q)
B为长pq的列矩阵
function B=rightmatrix(x,p,y,q,z)
这样就可以算出列矩阵a_n,然后转化成aa

具体操作:
先输入x=…y=…z=…p=…q= (注意x,y,z是行向量)
输入要求的坐标 xx=? yy=?
然后输入
>> main(x,p,y,q,z,xx,yy)
如果把main 改成一般的程序,而不是函数形式,运行后就可以在Workspace里看到aa系数矩阵,她可能对你非常重要。
问题:权函数其实在拟合的函数非常重要,不过她只能按你遇到的具体问题来取,我这里为1。
      当p,q 越大时,拟合的函数与原数据的方差越小,但是有可能函数本身抖动非常厉害,可以做个图看出来。

Matlab中矩阵求逆(A/B或者inv(A)*B)会出现很多warning,不知道是否对结果有影响。

panlei 发表于 2009-6-27 12:28

程序原理

程序原理见附件,谢谢你了!

波波球 发表于 2009-6-28 17:11

不妨参考一下此贴
http://forum.vibunion.com/forum/viewthread.php?tid=1005

panlei 发表于 2009-7-1 08:26

哪位大侠帮帮忙!跪谢!!!!!!!!
页: [1]
查看完整版本: 求助-曲面拟合