2023-06-17《计算方法》- 陈丽娟 - 插值法（二）.md

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

Matlab计算方法插值法埃尔米特插值分段低次插值三次样条插值

一、埃尔米特插值

埃尔米特插值即是找到一个插值函数, 使得不仅在节点上与原函数值相同，且还要求与原函数在节点处有相同的一阶(n阶)导数。下面是该问题的数学表达：
对于节点, 寻找插值多项式使得

为了求解上述函数，可设两组函数分别满足
(1) 至多都是次多项式
(2) 以下条件成立：

以及

然后可设：

由于在处函数值与导数值均为0，因此必定含有, 则令, 其中为拉格朗日基函数，定义为：

根据的要求，立即可得到的取值：

即

同样的可以得到

最终得到的埃尔米特插值公式即是

其中

注意到的形式，对于求导运算来说是比较复杂难计算的，因此埃尔米特插值公式常用于较少次数的插值。
最后，我们给出一个埃尔米特插值的程序，其中对导数的处理则是用的近似值:

1. %Hermite插值法，传入待求x，已知节点xi,yi,yi_diff, 导数参数delta 
   

   ---

2. function [y] = Hermite(x, xi, yi, yi_diff, delta) 
   

   ---

3. y = 0; %初始化y 
   

   ---

4. n = length(xi); 
   

   ---

5. for i = 1:n 
   

   ---

6. base_insert = BaseInsert(x, i, n, xi); 
   

   ---

7. base_insert_diff = (BaseInsert(x+delta, i, n, xi) - BaseInsert(x, i, n, xi)) / delta; 
   

   ---

8. y = y + (1 - 2 * (x - xi(i)) .* base_insert_diff) .* base_insert .^2 ... 
   

   ---

9. * yi(i) + (x - xi(i)) .* base_insert .^2 * yi_diff(i); 
   

   ---

10. end 
    

    ---

11. end 
    

    ---

12.  
    

    ---

13. function z = BaseInsert(x, i, n, xi) 
    

    ---

14. z = 1; 
    

    ---

15. for j = 1:n 
    

    ---

16. if j ~= i 
    

    ---

17. z = z .* (x - xi(j)) ./ (xi(i) - xi(j)); 
    

    ---

18. end 
    

    ---

19. end 
    

    ---

20. end 
    

    ---

测试程序

1. %Hermite插值法，传入待求x，已知节点xi,yi,yi_diff, 导数参数delta 
   

   ---

2. xi = 1:5; 
   

   ---

3. yi = sin(xi); 
   

   ---

4. yi_diff = cos(xi); 
   

   ---

5. delta = 1e-4; 
   

   ---

6. x = 1:0.1:5; 
   

   ---

7. y = Hermite(x, xi, yi, yi_diff, delta); 
   

   ---

8. plot(x,y) 
   

   ---

9. hold on 
   

   ---

10. scatter(xi,yi) 
    

    ---

给定五个点的插值结果：

给定10个点的插值结果：

可以看到由于埃尔米特插值的次数为, 因此相比于拉格朗日和牛顿插值法所得到的曲线不太理想。

二、分段低次插值

为了解决上述提到的高次插值结果不理想问题，考虑每次只插值少数节点，得到分段低次插值法。

由于分段线性插值在节点处不平滑，下面主要考虑分段埃尔米特插值。
设节点上的函数值为, 分段埃尔米特插值函数 满足

1. 为区间上一阶导数连续的函数。
2. 在各节点处函数值和其导数均与原函数相等。
3. 在每个区间上面是三次多项式。

按照埃尔米特插值公式可以直接给出分段的结果：

最后插值结果为

上述插值方法在左右两个端点处会强制变为0，此处不知道是否是理解错误
由上式我们给出下述插值程序

1. %Hermite插值法，传入待求x，已知节点xi,yi,yi_diff, 
   

   ---

2. function [y] = SeqHermite(x, xi, yi, yi_diff) 
   

   ---

3. y = 0; %初始化y 
   

   ---

4. n = length(x); 
   

   ---

5. for i = 1:n 
   

   ---

6. %先找到x对应的位置k 
   

   ---

7. for k = 1:n 
   

   ---

8. if x(i) >= xi(k) && x(i) <= xi(k+1) 
   

   ---

9. break 
   

   ---

10. end 
    

    ---

11. end 
    

    ---

12. y(i) = yi(k) * Calh(x(i), k, xi) + ... 
    

    ---

13. yi(k+1) * Calh(x(i), k+1, xi) + ... 
    

    ---

14. yi_diff(k) * CalH(x(i), k, xi) + ... 
    

    ---

15. yi_diff(k+1) * CalH(x(i), k+1, xi); 
    

    ---

16. % y(i) = 0; 
    

    ---

17. % for j = 1:length(xi) 
    

    ---

18. % y(i) = y(i) + yi(j) * Calh(x(i), j, xi) + yi_diff(j) * CalH(x(i), j, xi); 
    

    ---

19. % end 
    

    ---

20. end 
    

    ---

21. end 
    

    ---

22.  
    

    ---

23. function [hix] = Calh(x, k, xi) 
    

    ---

24. if k == 1 || k == length(xi) 
    

    ---

25. hix = 0; 
    

    ---

26. else 
    

    ---

27. x1 = xi(k-1); 
    

    ---

28. x2 = xi(k); 
    

    ---

29. x3 = xi(k+1); 
    

    ---

30. if x >= x1 && x <= x2 
    

    ---

31. hix = ((x-x1)/(x2-x1))^2 * (1 + 2 * (x - x2)/ (x1 - x2)); 
    

    ---

32. elseif x >= x2 && x <= x3 
    

    ---

33. hix = ((x-x3)/(x2-x3))^2 * (1 + 2 * (x - x2)/ (x3 - x2)); 
    

    ---

34. else 
    

    ---

35. hix = 0; 
    

    ---

36. end 
    

    ---

37. end 
    

    ---

38.  
    

    ---

39. end 
    

    ---

40.  
    

    ---

41.  
    

    ---

42. function [Hix] = CalH(x, k, xi) 
    

    ---

43. if k == 1 || k == length(xi) 
    

    ---

44. Hix = 0; 
    

    ---

45. else 
    

    ---

46. x1 = xi(k-1); 
    

    ---

47. x2 = xi(k); 
    

    ---

48. x3 = xi(k+1); 
    

    ---

49. if x >= x1 && x <= x2 
    

    ---

50. Hix = ((x-x1)/(x2-x1))^2 * (x - x2); 
    

    ---

51. elseif x >= x2 && x <= x3 
    

    ---

52. Hix = ((x-x3)/(x2-x3))^2 * (x - x2); 
    

    ---

53. else 
    

    ---

54. Hix = 0; 
    

    ---

55. end 
    

    ---

56. end 
    

    ---

57. end 
    

    ---

测试程序

1. %Hermite插值法，传入待求x，已知节点xi,yi,yi_diff, 导数参数delta 
   

   ---

2. xi = 1:10; 
   

   ---

3. yi = sin(xi); 
   

   ---

4. yi_diff = cos(xi); 
   

   ---

5. delta = 1e-4; 
   

   ---

6. x = 1:0.1:10; 
   

   ---

7. y = Hermite(x, xi, yi, yi_diff, delta); 
   

   ---

8. plot(x,y) 
   

   ---

9. hold on 
   

   ---

10. scatter(xi,yi) 
    

    ---

11.  
    

    ---

12. y = SeqHermite(x, xi, yi, yi_diff); 
    

    ---

13. plot(x,y) 
    

    ---

14. hold on 
    

    ---

15. scatter(xi,yi) 
    

    ---

最终得到下面的结果

可以明显看到在左右端点处没有和原始值相同，但是在内部的结果则比较理想。

三、Matlab的插值命令

在Matlab中使用interp1()函数进行插值。
使用方法 , 其中

1. vq是返回的插值
2. x是原始节点x坐标
3. v是节点对应的函数值
4. xq是需要插值的点
5. method。在matlab的文档里面有method的详细介绍：
   

例子：

1. x = 1:10; 
   

   ---

2. v = cos(x); 
   

   ---

3. xq = 0:0.1:11; 
   

   ---

4. methods = ["linear", 'nearest', 'next', 'previous', 'pchip', 'cubic',... 
   

   ---

5. 'v5cubic', 'makima', 'spline']; 
   

   ---

6. m = length(methods); 
   

   ---

7. n = length(xq); 
   

   ---

8. vq = zeros(n,m); 
   

   ---

9. for i = 1:m 
   

   ---

10. vq(:,m) = interp1(x,v,xq,methods(i)); 
    

    ---

11. plot(xq, vq(:,m)); 
    

    ---

12. hold on 
    

    ---

13. end 
    

    ---

14. legend(methods) 
    

    ---

15. hold off 
    

    ---

matlab自带函数插值结果

四、两道程序设计习题的解答

10. 已知, , 用线性插值和三次样条插值求的值。

1. x = [0.1, 0.8, 1.3, 1.9,2.5,3.1]; 
   

   ---

2. y = [1.2, 1.6, 2.7, 2.0, 1.3, 0.5]; 
   

   ---

3. xq = interp1(x,y,2) 
   

   ---

4. xq = interp1(x,y,2,'spline') 
   

   ---

11. 1990年到2010年每隔10年的产量为75.995, 91.972, 105.711, 123.203, 131.699, 150.697, 179.323, 203.212, 226.505, 249.633, 256.344, 267.893, 给出每隔一年的三次样条插值曲线。

1. x = 1900:10:2010; 
   

   ---

2. y = [75.995, 91.972, 105.711, 123.203, 131.699, 150.697, 179.323, ... 
   

   ---

3. 203.212, 226.505, 249.633, 256.344, 267.893]; 
   

   ---

4. xq = 1900:2010; 
   

   ---

5. vq = interp1(x,y,xq,'spline') 
   

   ---

6. plot(xq,vq) 
   

   ---

7. hold on 
   

   ---

8. scatter(x,y) 
   

   ---

9. hold off 
   

   ---

题11