例题
给出的函数关系表，分别利用牛顿插值法计算
的近似值。
| 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | |
| 1.105171 | 1.221403 | 1.349859 | 1.491825 | 1.648721 |
参考代码一：Python代码实现（自编码）
import math
"""
:difference_quotient差商函数
"""
def difference_quotient(x=None):
if x is None:
x = list()
if len(x) < 2:
raise ValueError('X\'s length must be bigger than 2');
ans = 0
for i in range(len(x)):
temp = 1.0
for j in range(len(x)):
if j == i:
continue
temp *= (x[i]-x[j]);
ans += (f(x[i])/temp)
return ans
def Calculate(data_x, x):
ans = f(data_x[0])
print(ans)
if len(data_x) == 1:
return ans
else:
temp = 1;
for i in range(len(data_x)-1):
temp*=(x-data_x[i])
ans+=difference_quotient(data_x[:i+2])*temp
return ans
def f(x):
return math.exp(x)
x = [0.1, 0.2, 0.3, 0.4, 0.5];
y = [1.105171,1.221403,1.349859,1.491825,1.648721];
target_point = 0.27;
fitted = Calculate(x,target_point);
real_value = f(target_point);
xx = [target_point]+x;
f_n = difference_quotient(xx);
w = []
for n in range(len(x)):
w.append(target_point-x[n])
Remainder = f_n;
for n in range(len(x)):
Remainder*=w[n];
print(f"点{target_point}处的真实值为{real_value}");
print(f"点{target_point}处的插值结果为{fitted}");
print(f"牛顿插值余项为{Remainder}");Python编码计算结果
参考代码二：MATLAB代码（自编码）
%可运行代码
clc,clear
format long
x = [0.1,0.2,0.3,0.4,0.5];%可以根据题目更换
y = [1.105171,1.221403,1.349859,1.491825,1.648721];%可以根据题目更换
newpoint = 0.27;%可以根据题目更换
newtonfitted = CalculateNewTon(x,newpoint);
xx = [newpoint,x];
f_n = difference_quotient(xx(:));
w=[];
for num=1:length(xx)-1
w = [w,newpoint-x(num)];
end
Remainder = f_n;
for m = 1:length(x)
Remainder = Remainder*w(m);
end
disp('点0.27处的插值结果为');
disp(newtonfitted);
disp('点0.27处的牛顿插值余项为');
disp(Remainder);
format shortdifference_quotient函数
difference_quotient函数用于计算牛顿差商的函数
function result = difference_quotient(x)
if isempty(x)
x = [];
end
if length(x)<2
error('Error:length must be bigger than 2');
end
ans0 = 0;
for i = 1:length(x)
temp = 1.0;
for j = 1:length(x)
if j==i
continue
end
temp = temp*(x(i)-x(j));
end
ans0 = ans0+(exp(x(i))/temp);
end
result = ans0;CalculateNewTon函数
CalculateNewTon函数用于计算最终插值结果
function result = CalculateNewTon(data_x,x)
ans1 = exp(data_x(1));
n = length(data_x);
if n == 1
result = ans1;
else
temp = 1;
for i = 1:length(data_x)-1
temp = temp*(x-data_x(i));
ans1 = ans1+difference_quotient(data_x(1:1:i+1))*temp;
end
result = ans1;
endMATLAB编码结果
附C++参考代码
#include<iostream>
#include<cmath>
#include<iomanip>
const int N = 10000;
using namespace std;
double f(double x)
{
return exp(x);
}
double difference_quotient(double x[],int count)
{
//count用来记录传入数组的大小
if(count<2){
throw "X\'s length must be bigger than 2";
}
double ans = 0;
for(int i=0;i<count;i++){
double temp=1.0;
for(int j=0;j<count;j++){
if(j==i) continue;
temp*=(x[i]-x[j]);
}
ans+=(f(x[i])/temp);
}
return ans;
}
double CalculateNewton(double data_x[],double x)
{
int count=0;
double second_data[N];
double result = f(data_x[0]);
double temp = 1.0;
for(int i=0;i<4;i++)
{
temp*=(x-data_x[i]);
for(int k=0;k<i+2;k++) {
second_data[k]=data_x[k];
count=k;
}
result+=difference_quotient(second_data,count+1)*temp;
}
return result;
}
int main()
{
double x[] = {0.1,0.2,0.3,0.4,0.5};
double y[] = {1.105171,1.221403,1.349859,1.491825,1.648721};
double point = 0.27;
double fitted = CalculateNewton(x,point);
double real_value = f(point);
double xx[]={0.27,0.1,0.2,0.3,0.4,0.5};
double f_n = difference_quotient(xx,6);
double w[N];
for(int i=0;i<5;i++){
w[i]=(point-x[i]);
}
double remainder = f_n;
for(int i=0;i<5;i++)
{
remainder*=w[i];
}
cout<<"真实值为"<<fixed<<setprecision(15)<<real_value<<endl;
cout<<"插值结果为"<<fixed<<setprecision(15)<<fitted<<endl;
cout<<"拉格朗日插值余项"<<fixed<<setprecision(15)<<remainder<<endl;
return 0;
}C++编码执行结果
参考书目
[1] 李新栋，许文文，张绪浩，任永强. 基于Python的计算方法[M]. 北京：电子工业出版社，2023.
