曲面拟合（可以得到函数）

glise

对 n 个三维坐标（x,y,z) ， 拟合成这样一个函数

Z(x,y)= aa(i,j)*x^i*y^j,   i 从0-->p ,j 从0-->q 求和。

具体的看里面的readme，说的很详细

1. function A=leftmatrix(x,p,y,q)
   
2. % A*a=B  a 即为系数列矩阵  
   
3. % A为左边(p-1)(q-1) 乘 (p-1)(q-1) 的矩阵
   
4. % x,y 为长度一样的列矩阵 也就是给定离散点的x,y坐标
   
5. % p,q为拟合的函数中x,y的指数
   
6. 
   
7. m=length(x);
   
8. if (nargin~=4) & (m~=length(y)), error('error check check!'); end
   
9. 
   
10. A_length=p*q;                        % A 为p*q阶的方阵
    
11. A=zeros(A_length,A_length);          % 赋值0
    
12. for i=1 : p*q
    
13.     for j= 1 : p*q
    
14.         x_z=quotient(j-1,q)+quotient(i-1,q);     % x 的指数   quotient为求商
    
15.         y_z=mod(j-1,q)+mod(i-1,q);            % y 的指数
    
16.         A(i,j)=qiuhe(x,x_z,y,y_z);            
    
17.     end
    
18. end

复制代码

1. function he=qiuhe(x,p,y,q,z)
   
2. % he    x^p*y^q 从1－>m的和
   
3. % x,y 行向量 长度相同; p,q 为x,y系数
   
4. %x=[2 3]; y=[ 3 4]; p=2; q= 2;
   
5. 
   
6. m=length(x);
   
7. if (nargin<4 )&(m~=length(y))            %输入量至少为四，x,y行向量长度必需一样
   
8.     error('error check check!');
   
9. end
   
10. if nargin==4, z=ones(m,1); end  %没有 z , 默认为单位行向量
    
11. 
    
12. he=0;
    
13. for i=1:m
    
14.     he=he+x(i)^p * y(i)^q*z(i);     % 1-->m 求和
    
15. end

复制代码

1. function sh=quotient(x,y)
   
2. % sh 为 x/y 的商
   
3. sh=(x-mod(x,y))/y;

复制代码

1. function B=rightmatrix(x,p,y,q,z)
   
2. % A*a=B
   
3. % B  为一个列向量 长为p*q
   
4. % x y z 为点的坐标 , p q 为x y指数
   
5. 
   
6. if nargin~=5,  error('error check check! rightmatrix'); end
   
7. B=zeros(p*q,1);
   
8. for i=1 : p*q
   
9.     x_z=quotient(i-1,q); y_z=mod(i-1,q); B(i,1)=qiuhe(x,x_z,y,y_z,z);
   
10. end

复制代码

1. function ff=main(x,p,y,q,z,xx,yy)
   
2. % x y z 坐标向量 长度要一样
   
3. % p ,q 为拟合函数中x,y 的系数
   
4. % xx yy 为 需要拟合的数据 给出(x,y) 坐标 求z
   
5. 
   
6. A=leftmatrix (x,p,y,q);             % A*a_n=B
   
7. B=rightmatrix(x,p,y,q,z);
   
8. 
   
9. %a_n=inv(A)*B;
   
10. a_n=A\B;                          % 求a_n    inv(A)*B 效果不好 （存疑）
    
11. 
    
12. for i=1 : length(a_n)             % 把长为p*q 的 a_n 列向量 转还 成p x q  的 aa 矩阵
    
13.     ii=quotient(i-1,q)+1;         % quotient求商
    
14.     jj=mod(i-1,q)+1; aa(ii,jj)=a_n(i,1);
    
15. end
    
16.                                  
    
17. ff=0;                              % ff 是 xx,yy 带入所拟合的函数 求出 z
    
18. for i=1 : p                        % 函数为 aa(i,j)*x^i*y^j  (i=0...p  ,j=0...q)
    
19.     for j=1 : q                    % aa 为系数  p x q 的矩阵
    
20.         ff=ff+aa(i,j) * xx^(i-1) * yy^(j-1);
    
21.     end
    
22. end

复制代码

[ 本帖最后由 ChaChing 于 2009-4-7 13:12 编辑 ]

goldenjay

为什么这个程序放到matlab里面会有问题？<BR>

galary

有没有对三维自变量情况下的拟合或者回归之类的方法？
我想请问，如果我的数据是一个四维的呢，就是测量数据受到三个因素的影响，即d=f(a,b,c),能拟合出函数来么？用什么方法解决？

[ 本帖最后由 ChaChing 于 2010-4-7 20:25 编辑 ]

chbo76

好，不错。
谢谢。支持一下。:lol

sunbader

x,y长度必须一样，这也太苛刻了！

lihuafeng

程序没任何问题，帮了我大忙，非常感谢！

guan_2984

xx,yy应该分别输入什么啊?

2365215

收藏了，谢谢lz

luxyer

这个程序具体该怎么运行呢？

wpaul

先感谢一下！谢谢。

ChaChing

原帖由 sunbader 于 2008-5-26 13:49 发表
x,y长度必须一样，这也太苛刻了！

x,y长度可以不一样吗?

lamb

这要不能用，我就真鸟了

misslfx

% xx yy  为需要拟合的数据 给出(x,y) 坐标 求z

xx,yy可以是你想要插值加密的点吧

[ 本帖最后由 ChaChing 于 2010-4-7 20:30 编辑 ]

ChaChing

怀疑没看程序说明就提问!? LS真体贴好心
还有三人问过(后面两人的便删除了!)

cfdbear

很好，感谢！！！！