0、回顾:对于 拉格朗日插值多项式与牛顿插值多项式的统一
在插值节点、插值条件相同的情况下,二者本质一样,只是计算过程不一样。
牛顿插值适合需要增加节点,提高精度的情况,不需要重新开始计算,可以利用上一步的计算结果。而可以利用上一步的计算结果的特性也是牛顿插值比拉格朗日插值较好的一点。即牛顿插值具有承继性。
问题:是不是节点越多,结果精度越高呢? 肯定不是这样的。
举个例子:函数
1、将 5等分,以分点做插值节点。
为了图像化,所以使用matlab来实现牛顿法插值。代码附在最后。
2、将 10等分,以分点做插值节点。
可以看出,当n增大时,图像在[-1,1]附近有着激烈振荡现象,这种现象称为龙格现象。
matlab代码:
调用函数
function[y,A,C,L] = newtonInter(X,Y,x)
n = length(X);
m = length(x);
for t = 1 : m
z = x(t);
A = zeros(n,n);
A(:,1) = Y';
for j = 2 : n
for i = j : n
A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1)); %差商矩阵
end
end
C = A(n, n);
for k = (n-1):-1:1
C = conv(C, poly(X(k)));
d = length(C);
C(d) = C(d) + A(k,k);
end
y(t) = polyval(C,z);
end
L = poly2sym(C);
执行文件
clear,clc;
% X = [-1, -0.5, 0, 0.5, 1];
% Y = [0.0385, 0.1379, 1, 0.1379, 0.0385];
X = [-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1];
Y = [0.0385, 0.0664, 0.1379, 0.39, 1, 0.39, 0.1379, 0.0664, 0.0385];
x = linspace(-1,1,250);
[y,A,C,L] = newtonInter(X, Y, x);
y1 = 1 ./ (1 + 25 * x.^2);
hold on
plot(X, Y, 'or', x, y, '.k',x,y1,'-b');
legend('样本点','牛顿插值估算','1/(1+25x^2)');
1、分段线性插值
对于龙格现象,计算是在整个区间上的,如果想解决这个问题,那么就把整个区间分成一段一段的来计算,这就是分段线性插值。
设 n+1 个节点
在区间上做线性插值,,在图形上即为把两点用线段相连,n条线段组成折线,该折线对应的函数称为分段线性插值函数。
插值基函数为:
首先要确定间隔符号,使得
令 、 ,则
插值结果为:(节点越多,越准确)
可以看出,分段线性插值也是存在缺点的,那就是连线是折线,整体的曲线不够光滑。
matlab代码:
function v = piecelin(x,y,u)
delta = diff(y) ./ diff(x);
n = length(x);
k = ones(size(u));
for j = 2:n-1
k( x(j) <= u) = j;
end
s = u - x(k);
v = y(k) + s.*delta(k);
end
执行代码;
clear,clc;
X = [-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1];
Y = [0.0385, 0.0664, 0.1379, 0.39, 1, 0.39, 0.1379, 0.0664, 0.0385];
x = linspace(-1,1,250);
y = piecelin(X,Y,x);
plot(X, Y, 'or', x, y, '.k');
legend('样本点','分段线性插值估算');
C++代码 (其计算出的点与matlab一致):
#include<iostream>
#include<vector>
#include<fstream>
using namespace std;
bool piecelin(double X[], double Y[], vector<double> &x, vector<double> &v);
int main()
{
double x[] = { -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1 };
double y[] = { 0.0385, 0.0664, 0.1379, 0.39, 1, 0.39, 0.1379, 0.0664, 0.0385 };
vector<double> u;
for (double i = -1; i <= 1; i = i + 0.01) {
u.push_back(i);
}
vector<double> v;
/*piecelin(x, y, u,v);*/
int n = sizeof(x) / sizeof(x[0]);
int k = u.size();
for (int i = 0; i < n - 1; i++) {
for (int j = 0; j < k; j++) {
if (x[i] <= u[j] && u[j] < x[i + 1]) {
double s = u[j] - x[i];
double delta = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
v.push_back(y[i] + s * delta);
}
if (u[j] > x[i + 1])
break;
}
}
int size = v.size();
std::ofstream output;
output.open("piecelin.txt");
for (unsigned i = 0; i < size; i++) {
output << u[i] << '\t' << v[i] << std::endl;
}
output.close();
cout << "Hello World !" << endl;
}