欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

一元以及二元多项式插值拟合(泰勒)

程序员文章站 2022-07-05 17:20:26
...

申明: 仅个人小记
根本上是基于泰勒公式,包括一元的和二元的泰勒定理。 泰勒用多项式逼近的思想。

效果展示

一元

一元以及二元多项式插值拟合(泰勒)

二元

一元以及二元多项式插值拟合(泰勒)

原理交代

一元

一元以及二元多项式插值拟合(泰勒)

二元

一元以及二元多项式插值拟合(泰勒)
其他推导部分和一元一样,本质上还是解线性方程组。

Matlab代码

一元
% 本质上就是n个方程解n个未知数,这里的未知数是待求函数的所有系数
% Ac=Y   A是由X组成的范德蒙德行列式,根据范德蒙德行列式的性质,
%为保证可解,X中不允许出重复的数值

X = 1:10;
Y = [4,5,1,8,2,-1,6,7,4,11];

% 很有意思,用到了范德蒙德行列式
% 因为范德蒙德行列式有很好的技巧性的计算方法,所以可能提供更好的计算方法
% 因为是范德蒙德行列式,所以,很容易知道什么情况该行列式不为零

n = length(X);
A = ones(n,n);

for j = 2:n % 从第二列开始,根据X计算相应的范德蒙德矩阵
   for i = 1:n
      A(i,j) = X(i)*A(i,j-1); 
   end
end

c = inv(A)*Y'; % 得到系数c, 即得到了相应的拟合函数
% f(x) = c0+c1*x+c2*x^2+...+cn-1 * x^(n-1)
% 下面绘出拟合函数
x = min(X):0.1:max(X); %
y = zeros(1,length(x));

for i = 1:length(x) % 带入x 计算 y
    t = 1;
    for j = 1:n
        y(i) = y(i)+t*c(j);
        t = t * x(i);
    end
end

subplot(311)
plot(X,Y,'r')
title('原数据点')
subplot(312)
plot(x,y,'g')
hold on
plot(X,Y,'O')
title('拉格朗日插值结果')
hold off
subplot(313)
plot(x,y,'g')
hold on
plot(X,Y,'r')
plot(X,Y,'ro')
title('数据比对')
hold off
二元
close all
clear all
X = [1 2 3 4 5 6 7 8 9 10]
Y = [6 2 3 12 9 9 7 3 1 9];
Z = [3 2 5 6 3 9 11 9 8 12];

n = length(X); % 必须保证 n = 1+2+3+...+m, m为整数
m = floor(sqrt(2*n))-1; % 计算相应目标函数的阶数, 从 0 阶开始

%% 数据计算准备
% tt 中的内容及意义
%         0 阶   1阶     2阶      3阶   
% x的次幂  0   , 0 1 ,  0 1 2 ,  0 1 2 3 , ...
% y的次幂  0   , 1 0 .  2 1 0 ,  3 2 1 0 , ...
tt = zeros(2,n);
k = 1;
for i = 0:m
    for j = 0:i
        tt(1,k) = j;
        tt(2,k) = i-j;
        k = k+1;
    end
end

%% 根据tt, X, Y, 计算相应的系数矩阵 A
A = ones(m,m);
for i = 1:n
    k = 1;
    for j = 1:n      
        A(i,j) = power(X(i),tt(1,k))*power(Y(i),tt(2,k));
        k = k+1;
    end
end

c = inv(A)*Z'; % 得到目标函数的系数, 即得到 z = f(x,y)

%% 绘制目标拟合函数图
% z = f(x,y)

[x, y] = meshgrid(min(X):0.5:max(X),min(Y):0.5:max(Y));

% 计算z值
z = zeros(size(x,1),size(x,2)); % 只是赋予z和x同样的规格
for i = 1:size(x,1)
    for j = 1:size(x,2)
        for k = 1:n
            z(i,j) = z(i,j) + c(k)*power(x(i,j),tt(1,k))*power(y(i,j),tt(2,k));
        end
    end
end
subplot(211)
mesh(x,y,z)
xlabel('X')
ylabel('Y')
zlabel('Z')
hold on
plot3(X,Y,Z,'ro')
hold off

subplot(212)
mesh(x,y,z)
xlabel('X')
ylabel('Y')
zlabel('Z')
hold on
plot3(X,Y,Z,'ro')
plot3(X,Y,Z,'r')

2018年1月24日 13:41:16 Written by Jack Lu