课本上给出的那种方法请参考帖子“(Newton插值多项式——cos x)帮忙看下这程序为啥不对好吗?谢谢老师和同学”,不过下面的程序在计算差商或差分时更节省存储空间。
1 cos(x) (未给出误差分析,请参考插值多项式的误差余项个人完成)
2 Bessel函数
1 cos(x) (未给出误差分析,请参考插值多项式的误差余项个人完成)
- 代码:
% P37 例3 函数cos(x)的Newton插值逼近
clear; clc; close all;
format long;
X=[0.0 0.2 0.4 0.6 0.8 1.0 1.2];
Y=[1.0 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358];
n = length(X);
A = Y';
% 计算cos(x)的各阶差商,存入A
% 在X(1),X(2),..,X(k)处的k-1阶差商存入A(k), k=1,2,...,n
for i = 1:n-1
for j = n:-1:i+1
A(j) = ( A(j) - A(j-1))/(X(j) - X(j-i));
end
end
A
% 计算Newton插值多项式在x0处的函数值y0
x0 = 0.3;
y0 = A(1);
temp = 1;
for j = 2:n
temp = temp * (x0 - X(j-1));
y0 = y0 + A(j)*temp;
end
% 输出
x0
y0
2 Bessel函数
- 代码:
% P44 例4 Bessel函数的Newton向前插值逼近
clear;clc;
format long
x = 2.4:0.1:2.9;
y = besselj(0,x);
n = length(x);
% 计算各阶差分。通过下面的计算,将x0=2.4处的第i阶差分存入y(i+1)
for i = 1:n-1
for j = n:-1:i+1
y(j) = y(j)-y(j-1);
end
end
% 通过Newton向前插值公式计算2.45处的函数值
f = y(1);
t = 0.5; % 0<t<1, 步长h = 0.1, t=0.5即现在考虑的点为2.45 + t*h = 2.45
co = 1;
for i = 1:n-1 % i 表示差分阶数
% co 表示组合数(t,i)
co = co * (t-i+1)/i;
% 输出 利用从点2.4及后面i个节点得到插值多项式 在2.45处的函数值
f = f + co*y(i+1)
end
% 2.45处精确函数值
f_exact = besselj(0,2.45)