|
申明： 仅个人小记 根本上是基于泰勒公式，包括一元的和二元的泰勒定理。 泰勒用多项式逼近的思想。
效果展示
一元
二元
原理交代
一元
二元
 其他推导部分和一元一样，本质上还是解线性方程组。
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
|