想跟大家请教一个我遇到曲线拟合的问题

zhiyucn

我这里有3组数据x1, x2 , y 三组数据都是 2000行1列的结构，我的任务是找出y与x1,x2 之间的关系, 就是 y = f(x1,x2)。

我尝试过用多项式，可是一直没找到合适的。现在我在尝试用rational Funktion

比如说   y_fit = (a1*x1+constant) / (b2*x2^2+b1*x2-1)               

C = [b2.*x2.^2.*y +b1.*x2.*y  –a1.*x1  -ones(size(x2)) ];     %求系数

W = C\y;

Y = [b2.*x2.^2.*y +b1.*x2.*y  –a1.*x1  -ones(size(x2)) ]*W;  %验证是否正确

我求出了相应的系数后，经过了验证，结果Y和y的结果吻合的很好。 然后我又将系数w带回上面的函数中，想再验证一下结果是否正确，可是出来了完全不同的结果。请教各位不知道你们有没有遇到这种情况，请问是什么地方出现错误了呢? 先谢谢了。

[ 本帖最后由 zhiyucn 于 2007-8-4 18:29 编辑 ]

aileute

为什么是C\y，而不是C./y呢？或者C.\y呢？

aileute

我有个问题，你为什么找不到合适的多项式呢？如果是内插值来重构模型的话，总是可以找到满足要求的多项式的呀！是不是你采用的pade拟合，对你来说物理上跟有意义些呢？因为你可以从中求出零极点，然后再反变换回时域什么的。
只是好奇，随便问问。

更正一下，仔细看了一下，你用的是又两个变量的有理分式，所以不是pade拟合。Sorry!

[ 本帖最后由 aileute 于 2007-8-4 12:44 编辑 ]

zhiyucn

aileute 已经给你发了短信，希望能够跟你电话交流，那样比较方便。

[ 本帖最后由 zhiyucn 于 2007-8-4 11:06 编辑 ]

zhiyucn

原帖由 aileute 于 2007-8-4 10:49 发表
我有个问题，你为什么找不到合适的多项式呢？如果是内插值来重构模型的话，总是可以找到满足要求的多项式的呀！是不是你采用的pade拟合，对你来说物理上跟有意义些呢？因为你可以从中求出零极点，然后再反变换回 ...

C\y 可以求出所要的系数矩阵，结构也是 2000 行1列的。如果用C./y 那就不是这样的结果了。

请教怎麽样用内插值重构模型?

aileute

我自己造了一组数据，做了一个测试，发现样本的数量在这个问题中起了很大作用。一下是测试程序
clc; clear; close all;
% Version: 04.08.2007
% Author: Wang,Lei
% This file is used to test a fitting problem
% y = (x1 + 1)/(x2^2 + 2*x2 + 1)
x1 = (0:0.01:10)'; % 我主要调整这里0.2，0.02，0.01，来调整样本数量
x2 = (0:0.01:10)';
y = (x1 + 1)./(x2.^2 + 2*x2 + 1);
A = [-y.*x2.^2  -y.*x2  x1 ones(size(x1))]; %你得这里的符号也写错了，可能也是得不到结果的原因。
% A = [y.*x2.^2 y.*x2  -x1 -1*ones(size(x1))];  % 你的系数矩阵
[U S V] = svd(A'*A);
Sinv = zeros(size(S));
for i = 1:size(S,1)
    Sinv(i,i) = 1/S(i,i);
end
W = inv(V')*Sinv*inv(U) * A'* y
y_fit = (W(3).*x1 + W(4))./(W(1).* x2.^2 + 2*x2 + 1);
% y_fit = A * W; % 你得测试程序

plot(1:length(y),y,'r*-', 1:length(y_fit), y_fit, 'bo-');

你运行一下以上这个程序就知道了。

zhiyucn

谢谢了你的程序，我会好好的研究一下的。

zhiyucn

我把3组数据重新调整了一下，让他们尽量接近我的数据。

clc; clear; close all;
% Version: 04.08.2007
% Author: Wang,Lei
% This file is used to test a fitting problem
% y = (x1*x2 + 1)/(x1*x2^2 + x1*x2 - 1)

x1 = (1e-10:1e-10:1e-8)';
x2 = 323*ones(size(x1));
y = (0.5:-0.008:-0.2992)' ;

A = [y.*x1.*x2.^2  y.*x1.*x2  -x1.*x2 -ones(size(x1))];

[U S V] = svd(A'*A);
Sinv = zeros(size(S));
for i = 1:size(S,1)
    Sinv(i,i) = 1/S(i,i);
end
W = inv(V')*Sinv*inv(U) * A'* y

y_fit = (W(3).*x1.*x2 + W(4))./(W(1).*x1.* x2.^2 +W(2).*x1.*x2 - 1);

semilogx(x1,y,'r*-', x1, y_fit, 'bo-');

结果如下

W =

  1.0e+005 *

   -0.0000
    0.0045
    2.4768
   -0.0000



[ 本帖最后由 zhiyucn 于 2007-8-4 20:01 编辑 ]