函数逼近——(Lagrange)拉格朗日插值法 | 北太天元 or Matlab

函数逼近——(Lagrange)拉格朗日插值法 | 北太天元 or Matlab

引用

CSDN

1.

https://blog.csdn.net/Math_Boy_W/article/details/137800068

拉格朗日插值法是一种常用的数值分析方法，用于通过已知数据点构造多项式函数，从而实现对未知数据点的预测。本文将详细介绍拉格朗日插值法的理论基础、算法实现以及在北太天元或Matlab中的具体应用。

一、Lagrange插值法

对于n + 1个样本点$(x_i, y_i), i=0,1,2,\dots ,n$，Lagrange插值多项式可以表示为：

$$
L_n(x)=\sum_{i=0}^{n}y_iI_i(x)
$$

其中，$I_i(x)$是Lagrange基函数，定义为：

$$
I_i(x)=\frac{(x-x_0)(x-x_1)\dots(x-x_{i-1})(x-x_{i+1})\dots(x-x_n)}{(x_i-x_0)(x_i-x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)}
$$

也可以表示为：

$$
I_i(x)=\frac{\omega_n(x)}{(x-x_i)\omega_n'(x_i)}
$$

其中，$\omega_n(x)=\prod_{j=0}^{n}(x-x_j)$。

$L_n(x)$是$f(x)$的n次多项式插值的Lagrange公式，也称为Lagrange插值多项式。

二、算法

♡Lagrange插值法的实现步骤如下：

1. 输入已知插值点的横坐标向量$x_0$和纵坐标向量$y_0$
2. 计算每个$x$对应的Lagrange基函数$I_i(x)$
3. 输入所求点的横坐标向量$x$
4. 结合$I_i(x)$与$y_0$计算每个$x$对应的$L_n(x)$
5. 输出对应的$L_n$

三、北太天元 or Matlab 实现

以下是Lagrange插值法的Matlab实现代码：

function  [y] = Lag_interp_v1(x0,y0,x)
    n1 = length(x0);   % n1表示样本点的个数 
    I = zeros(1,n1);   % 预留出要用的空间
    n2 = length(x);    % n2表示所求点的个数
    Ln = zeros(1,n2);
    for j=1:1:n2     %  依次代入 自变量 x(j) 
        omega_x = x(j)-x0;    % 数 - 矩阵 ，表示 [x(j) - x0(1),x(j)-x0(2),...,x(j)-x0(n1)]
        for i = 1:1:n1   %  对于x(j)求对应的基函数 I（i）
            w = x0(i) - x0;  % 同样是 数 - 矩阵
            % 这里使用 if 和 内置的 prod 函数代替了 for 循环
            % prod 表示矩阵内所有元素的乘积
            if i == 1    
                % omega_x(i+1:n1)表示向量的节选，第i+1个到第n1个元素
                I(i) = prod(omega_x(i+1:n1))/prod(w(i+1:n1));
            elseif i == n1
                I(i) = prod(omega_x(1:i-1))/prod(w(1:i-1));
            else
                I(i) = prod(omega_x(1:i-1))/prod(w(1:i-1));
                I(i) = I(i) * prod(omega_x(i+1:n1))/prod(w(i+1:n1));
            end 
        end
        %使用矩阵的乘积，行向量 × 列向量 得到一个值
        Ln(j) = y0*I'; 
    end    
    y = Ln;
end

将上述代码保存为Lag_interp_v1.m文件。

四、数值算例

利用$f(x) = \ln x$的如下数据：

进行Lagrange插值:

1. 计算$x = [0.412, 0.511, 0.666]$处的近似值
2. 计算$x_i = 0.3 + ih, h=0.01, i = 0, 1, 2, \dots, 50$处的近似值，并作图

调用函数[y] = Lag_interp_v1(x0,y0,x)，相应的实现代码为：

clc, clear all, format long;
x0 = linspace(0.4, 0.7, 4); % 输入样本点的横坐标
y0 = [-0.916291, -0.693147, -0.510826, -0.35765]; % 输入样本点的纵坐标
% 简单的算几个点
x1 = [0.412, 0.511, 0.666];
y1 = Lag_interp_v1(x0, y0, x1);
for i = 1:length(x1)
    fprintf('f(%f) = %f \n', x1(i), y1(i));
end
% 利用很多的点来画图
x2 = linspace(0.3, 0.8, 51); % 共51个要求的点
% 利用写好的 Lag_interp_v1 函数计算要求点的纵坐标
y2 = Lag_interp_v1(x0, y0, x2);
delta = abs(y2 - log(x2));
%作图
figure(1); %画出第一个图像
plot(x2, y2, 'b');
figure(2);
plot(x2, y2, 'b');
hold on
plot(x2, log(x2), 'r'); % y = lnx 的图像
hold off
% Ln(x) 与 lnx 的 误差
figure(3); % 画出第二个图像
plot(x2, delta, 'g');
% 文字形式表示出来
for j = 1:51
    fprintf('f(%f) = %f \n', x2(j), y2(j));
end

将上述代码保存为LagTest.m，运行后得到以下结果：

f(0.412000) = -0.886972 
f(0.511000) = -0.671305 
f(0.666000) = -0.407185 
f(0.300000) = -1.191936 
f(0.310000) = -1.161676 
f(0.320000) = -1.132046 
f(0.330000) = -1.103035 
f(0.340000) = -1.074630 

此处省略
f(0.770000) = -0.261515 
f(0.780000) = -0.248246 
f(0.790000) = -0.235059 
f(0.800000) = -0.221941 

$L_n(x)$的图像如下：

$L_n(x)$与$\ln(x)$的图像如下：

$L_n(x)$与$\ln(x)$的误差如下：