matlab：Lagrange插值函数构造

$Lagrange$基函数及$Lagrange$多项式的构造

$构造Lagrange基函数$

这里我们需要注意到一个小技巧，即在$Matlab$中，求和函数$sum$与求积函数$prod$先按列计算的，我们在计算行向量与列向量时体现不明显，但我们计算矩阵时，返回的结果是行向量，也就是每一列元素求和（求积）的结果构成的行向量，如下：

矩阵	求和	求积

所以我们对矩阵求行和时，需要将矩阵转置然后再利用$sum$函数。

在$Lagrange$基函数的构造中，我们知道$L_j(x)={\prod }_{i=0,i\ne j}^n\frac{(x-x_i)}{(x_j-x_i)}$。

令

$A=x-\left[\begin{array}{c}{x}_1\\ x_2\\ \dots \\ x_{i-1}\\ x_{i+1}\\ \dots \\ x_n\end{array}\right],B=\left[\begin{array}{c}{x}_j-x_1\\ x_j-x_2\\ \dots \\ x_j-x_{i-1}\\ x_j-x_{i+1}\\ \dots \\ x_j-x_n\end{array}\right]$

从而，$L_j(x)=\frac{\prod (A)}{\prod (B)}$，其实这里面，不论$B$是行向量还是列向量，都可以计算出正确的结果，但我们利用求积函数构造关于$x$的函数时，则会出现下列这几种情况；这与$matlab$中$prod$函数的构造原理是分不开的。这里为了看清我们简单的举个例子：

	行向量	列向量
	$x-[1,2,3]\iff [x-1,x-2,x-3]$	$x-[1;2;3]\iff [x-1;x-2;x-3]$

	行向量		列向量
求积	类型1：x输入为行向量	类型2：x输入为列向量	类型3：x输入为行向量	类型4：x输入为列向量

这里面不难发现，$Matlab$的$prod$函数是以按列求积的方式进行计算的。

1. 对类型1，我们有$prod([4,5]-[1,2,3])$时，做减法时向量维度不匹配，故报错；
2. 对类型2，我们有$prod([4;5]-[1,2,3])=prod(\left[\begin{array}{c}3,2,1\\ 4,3,2\end{array}\right])=[12,6,2]$;
3. 对类型3，我们有$prod([4，5]-[1;2;3])=prod(\left[\begin{array}{c}3,4\\ 2,3\\ 1,2\end{array}\right])=[6,24]$;
4. 对类型4，我们有$prod([4;5]-[1;2;3])$时，做减法时向量维度不匹配，故报错。

综上，只有类型2和类型4可以输出结果，但当输入变量为列向量类型2求积方式会改变为按列求积（这本身就是$prod$函数的运算方式），所以我们选择类型3，运算的结果符合$Language$差值基函数的要求（需要注意的是，我们的输入必须为行向量）。

function Lagrange_basis = Lagrange_basis(xi)
[a,b] = size(xi);
% 判断是否为矩阵，如果是则报错。
if a>1 && b>1
    error("Input error,please input vector")
end
% 判断xi是否为列向量，如果不是，则转置为列向量。
if a == 1
    xi = xi';
end
% Lagrange_basis将会以元胞数组的形式储存所有的Lagrange基函数
Lagrange_basis = cell(1,len);
for j = 1:len
    L_j = @(x)prod(x - xi([1:j-1,j+1:end]))./prod(xi(j)-xi([1:j-1,j+1:end]));
    Lagrange_basis{j} = L_j;
end

构造$Lagrange$函数

使用上述构造的基函数构造Lagrange函数

从而构造$Lagrange$函数是简单的，上面我们已经知道$Lagrange$基函数保存在一个元胞数组里面，接下来提取即可并于$y(i)$求积再求和即可；

function L = Lagrange_function(Lagrange_basis,y)
len = length(y);
% 初始化Language函数L=0
L = @(x) 0; 
% 累计求和
for i =1:len
    % 获取第i个Lagrange基函数
    Li = Lagrange_basis{i};
    % 第i个Lagrange基函数与y_i求积并与前面的部分Lagrange函数累加
    L = @(x) L(x) + Li(x) * y(i);
end

举例：

一个简单的例子，$y=x^4,x=[1,2,3,4,5],y=[1,16,81,256,625]$

函数图及Lagrange函数图	误差曲线

例子代码：

clc,clear
xi = [1:1:5]';
y =@(x) x.^4;
Lagrange_basis = Lagrange_basis(xi);
L = Lagrange_function(Lagrange_basis,y(xi));
x_v = 1:0.01:8;
plot(x_v,L(x_v),'k',x_v,y(x_v),'r')
legend('Lagrange插值函数','y=x^4')
figure
title('误差')
plot(x_v,y(x_v)-L(x_v))

• 需要注意的是，这里我们构造Lagrange基函数时需要使用的是列向量，而调用Lagrange函数(代码中的Lagrange_function函数)时需要使用行向量。

在构造基函数的同时构造加入Lagrange函数的构造

当然你如果不想存储Lagrange基函数，完全可以这样处理，我们沿用上面的例子，代码如下：

clc;clear
xi = [1:1:5]';
y =@(x) x.^4;
yi = y(xi);
Lagrange_fun = @(x) 0;
for i = 1:length(xi)
    Lagrange_fun = @(x)Lagrange_fun(x) + y(i)*(prod(x - xi([1:i-1,i+1:end]))./prod(xi(i)-xi([1:i-1,i+1:end])));
end
x_v = 1:0.01:8;
plot(x_v,Lagrange_fun(x_v),'k',x_v,y(x_v),'r')
legend('Lagrange插值函数','y=x^4')
figure
title('误差')
plot(x_v,y(x_v)-Lagrange_fun(x_v))

实验结果同上。

• 需要注意的是，这里我们构造Lagrange基函数时需要使用的是列向量，而调用Lagrange函数(代码中的Lagrange_function函数)时需要使用行向量。

利用矩阵或者向量的运算去掉for循环构造Lagrange函数

• 虽然我们确实构造出来了，但感觉这种写法并不好用。（只是做笔记的一部分而已）

• 去掉了Lagrange函数构造里面的for循环(用于基函数与对应y值求积后累加)，但是我们只是使用sum函数代替。

• 在这个代码中，需要注意的是，输入的插值节点必须为列向量，而对应输入的y值必须为行向量

• 最后，这个函数也就是Lagrangefun(x)一次只能求一个值，即x不能为向量，（这一点让人很难受），所以当我们的x为向量时，我们需要单独写一个循环。

举例：同上

实验结果也和上面一样。

函数图及Lagrange函数图	误差曲线

代码：

clc;clear
figure
xi = [1:1:5]';
y =@(x) x.^4;
yi = y(xi');
len = length(xi);
xi_m = xi .* ones(len,len); 
xi_m(1:len+1:end)=[];
xi_m = reshape(xi_m,[len-1,len]);
Lagrangefun = @(x) sum(yi.* prod(x'.*ones(1,len)-xi_m)./prod(xi'-xi_m));
x_v = 1:0.01:8;
y_v = zeros(1,length(x_v));
for i = 1:length(x_v)
   y_v(i) =  Lagrangefun(x_v(i));
end
plot(x_v,y_v,'k',x_v,y(x_v),'r')
legend('Lagrange插值函数','y=x^4')
figure
title('误差')
plot(x_v,y(x_v)-y_v)