2023-06-16《计算方法》- 陈丽娟 - 插值法（一）

合集 - 《计算方法》-陈丽娟(8)

1.2023-06-16 《计算方法》- 陈丽娟 - 绪论.md2023-06-16

2.2023-06-16《计算方法》- 陈丽娟 - 插值法（一）2023-06-16

3.2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法.md2023-06-274.2023-06-19《计算方法》- 陈丽娟 - 方程的近似解法（注解）2023-06-195.2023-06-18《计算方法》- 陈丽娟 - 方程的近似解法.md2023-06-186.2023-06-17《计算方法》- 陈丽娟 - 插值法（二）.md2023-06-177.2023-06-30《计算方法》- 陈丽娟 - 线性方程组的迭代解法.md2023-06-308.2023-06-28《计算方法》- 陈丽娟 - 向量和矩阵基础.md2023-06-28

收起

2023-06-16《计算方法》- 陈丽娟 - 插值法（一）

Matlab计算方法插值法

本章给出了一些基本的插值法理论和算法，附带解决部分习题。

一、拉格朗日插值

为了直观，这里部分符号和书中不一致，但是得到的形式更优美。

1. 一次拉格朗日插值

在给定区间上，已知, , 一次拉格朗日插值要求插值函数满足, .

由两点确定一条直线，显然可以得到如下的一个可行

改写为两点式，可得

上述函数可以看作是

其中称为线性插值基函数，且满足, , , 和.

2. 二次拉格朗日插值

根据上述一次拉格朗日插值方法，我们可以构建对于3个点和对应值的二次拉格朗日插值法。

首先给出二次拉格朗日插值线性插值基函数的条件：

由于三个点可以确定一条二次函数，设, 可以分别得到为

注意到, 带入上式可得

立即可得到二次拉格朗日插值公式

3. 多次拉格朗日插值

下面直接给出多次拉格朗日插值法的形式。
对于一个含有n个点和对应的值的拉格朗日插值问题，我们有下述解析函数

其中定义为(拉格朗日插值基函数)

唯一性定理

在次数不超过的几何多项式中，满足插值条件的插值多项式 是存在的，且是唯一的。

证明: 存在性由拉格朗日插值法可得。唯一性可假设存在是另一个满足条件的多项式，则对所有个点成立，但由次多项式最多只有个根得到矛盾。

最后在本节我们给出拉格朗日插值的Matlab程序：

1. % 拉格朗日插值，传入x_i, y_i, 和待求点x，返回x对应的插值和系数Lx 
   

   ---

2. function [y Lx] = LagInsert(x, xi, yi) 
   

   ---

3. %若x是数组，则对每个x返回值,以及系数组（如果需要返回的矩阵太大，则不返回） 
   

   ---

4. Lx = []; 
   

   ---

5. y = 0; 
   

   ---

6. n = length(xi); 
   

   ---

7. isReturnLx = true; 
   

   ---

8. if n * length(x) > 1e5 
   

   ---

9. isReturnLx = false; 
   

   ---

10. end 
    

    ---

11. for i = 1:n 
    

    ---

12. lx = 1; 
    

    ---

13. for j = 1:n 
    

    ---

14. if j ~= i 
    

    ---

15. lx = lx .* (x - xi(j))/ (xi(i) - xi(j)); 
    

    ---

16. end 
    

    ---

17. end 
    

    ---

18. y = y + lx * yi(i); 
    

    ---

19. if isReturnLx 
    

    ---

20. Lx = [Lx;lx]; 
    

    ---

21. end 
    

    ---

22. end 
    

    ---

23. end 
    

    ---

对进行测试

1. xi = 1:10; 
   

   ---

2. yi = sin(xi); 
   

   ---

3.  
   

   ---

4. x = 0:0.1:10; 
   

   ---

5. [y Lx] = LagInsert(x, xi, yi); 
   

   ---

6.  
   

   ---

7. plot(x,y) 
   

   ---

8. hold on 
   

   ---

9. scatter(xi,yi) 
   

   ---

得到结果

书中也有拉格朗日插值的代码，处理上有些微不同。

二、差商与牛顿插值

拉格朗日插值当插值点增加或者减少时，需要重新开始所有计算，牛顿法可以解决这个问题。

考虑对于个点的次多项式插值：

首先给出差商的定义：

差商

函数关于点, 的一阶差商; 关于的二阶差商; 以及阶差商

不知怎么的#F44336反正就有

1. 差商与节点的位置无关
   
   
2. 注意到位置无关性和差商定义，有
   
   
3. 差商与导数之间的关系
   
   
4. 若, 则.
5. 若, 则.

差商的计算-- 差商表。注意差商表中比如第四行第三列的计算方式为

我们给出差商表的计算程序：

1. %根据一组数据，计算对应的差商表(differience quotients) 
   

   ---

2. function [dq] = DQ(x, fx, dq) 
   

   ---

3. %x 是横坐标值 
   

   ---

4. %fx 是对应的函数值 
   

   ---

5. %首先判断dq是否是一张dq表 
   

   ---

6. if length(dq) == 1 
   

   ---

7. n = length(x); 
   

   ---

8. dq = zeros(n, n); 
   

   ---

9. dq(:,1) = fx; 
   

   ---

10. for i = 2:n %列 
    

    ---

11. for j = i:n %行 
    

    ---

12. dq(j,i) = (dq(j,i-1) - dq(j-1,i-1))/(x(j)-x(j-i+1)); 
    

    ---

13. end 
    

    ---

14. end 
    

    ---

15. else 
    

    ---

16. %在表后面添加值 
    

    ---

17. m = length(dq); 
    

    ---

18. n = length(x); 
    

    ---

19. dq_temp = dq; 
    

    ---

20. dq = zeros(n,n); 
    

    ---

21. dq(1:m,1:m) = dq_temp; 
    

    ---

22. dq(:,1) = fx; 
    

    ---

23. for i = 2:n 
    

    ---

24. for j = (max(m,i)):n 
    

    ---

25. dq(j,i) = (dq(j,i-1) - dq(j-1,i-1))/(x(j)-x(j-i+1)); 
    

    ---

26. end 
    

    ---

27. end 
    

    ---

28. end 
    

    ---

29. end 
    

    ---

该算法可以处理当差商表存在时快速计算新加入的点的差商。

得到差商表后，我们可以根据下述牛顿插值法：

结合差商程序，我们给出下列牛顿插值法的程序：

1. function [NI] = NewtonInsert(tx, x, fx, dq) 
   

   ---

2. %计算dq表 
   

   ---

3. [dq] = DQ(x, fx, dq); 
   

   ---

4. %计算插值的结果, 可以是向量 
   

   ---

5. NI = dq(1,1); 
   

   ---

6. for i = 2:length(dq) 
   

   ---

7. tq = 1; 
   

   ---

8. for j = 1:(i-1) 
   

   ---

9. tq = tq .* (tx-x(j)); 
   

   ---

10. end 
    

    ---

11. NI = NI + dq(i,i) * tq; 
    

    ---

12. end 
    

    ---

13. end 
    

    ---

最后依然给出实例代码：

1. % 测试生成的差商表 
   

   ---

2.  
   

   ---

3. x = 1:10; 
   

   ---

4. fx = sin(x); 
   

   ---

5.  
   

   ---

6. [dq] = DQ(x, fx, 1); 
   

   ---

7.  
   

   ---

8. % 增加x的数量 
   

   ---

9. y = [x x+0.1]; 
   

   ---

10. fy = sin(y); 
    

    ---

11.  
    

    ---

12. tx = （-1）:0.1:10; 
    

    ---

13. [NI1] = NewtonInsert(tx, x, fx, dq); 
    

    ---

14. [NI2] = NewtonInsert(tx, y, fy, dq); 
    

    ---

15. plot(tx,NI1) 
    

    ---

16. hold on  
    

    ---

17. plot(tx, NI2) 
    

    ---

18. hold on 
    

    ---

19. scatter(x, fx) 
    

    ---

插值结果如下图所示

可以看到在插值区间内的拟合效果很好，但是在插值区间外没有限制。