数值计算方法第4次作业

第四章

问题一

一、问题综述

在离地球表面高度为y处的重力加速度如下：

计算高度y=55000m处的重力加速度值。

二、问题分析

以高度y作为自变量，重力加速度的值为因变量。得到以下信息：

f(0)=9.8100;

f(30000)=9.7487;

f(60000)=9.6879;

f(90000)=9.6278;

f(120000)=9.5682;

本题要求的就是f(55000)的值。

以下将采用课堂中学到的Lagrange插值多项式法、Newton插值多项式法、分段低次插值法和样条插值法求解该问题。

三、问题解决

1. lagrange插值多项式法

对某个多项式函数，已知有给定的k+ 1个取值点：

其中对应着自变量的位置，而对应着函数在这个位置的取值。

假设任意两个不同的x j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：

其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为：

拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。

源程序lagrange.m

function [c,f]=lagrange(x,y,a)

% 输入：x是自变量的矩阵；y是因变量的矩阵；a是要计算的值的自变量；

% 输出：c是插值多项式系数矩阵；f是所求自变量对应的因变量；

m=length(x);

l=zeros(m,m); % l是权矩阵

f=0;

for i=1:m

v=1;

for j=1:m

if i~=j

v=conv(v,poly(x(j)))/(x(i)-x(j)); % v是l_i(x)的系数矩阵

end

end

l(i,:)=v; % l矩阵的每一行都是x从高次到低次的系数矩阵

end

c=vpa(y*l,10); % 对应阶次的系数相加，乘以y，显示10位有效数字

for k=1:m

f=f+c(k)*a^(m-k);

end

输入矩阵

x=[0 30000 60000 90000 120000]

y=[9.81 9.7487 9.6879 9.6278 9.5682]

a=55000

再运行源函数，可得：

c =

[ -2.057613169e-23, 4.938271605e-18, -3.703703702e-14, -0.000002046111111, 9.81]

f =

9.6979851723251649906109417384537

即此时得到的插值函数为：

f(x)=−2.06×10−23x4+4.94×10−18x3−3.70×10−14x2−2.05×10−6x+9.81所以可以得到本题的结果为：

f(55000)=9.6980

2. Newton插值多项式法

Newton插值多项式法通过各阶差商来得到插值函数。牛顿基本插值多项式为：N n(x)=f(x0)+f[x0,x1](x−x0)+⋯+f[x0,…,x1](x−x0)…(x−x n−1)其中

零阶差商 f[x i]=f(x i)

一阶差商 f[x i,x j]=f(x i)−f(x j) x i−x j

二阶差商 f[x i,x j,x k]=f[x i,x j]−f[x j,x k]

x i−x k

……

k阶差商 f[x i,x i+1,…,x i+k]=f[x i+1,…,x i+k]−f[x i,…,x i+k−1]

x i+k−x i

源程序newton.m:

function [c,d,f]=newton(x,y,a)

% 输入：x是自变量矩阵；y是因变量矩阵；a是要计算的值的自变量；

% 输出：c是插值多项式系数矩阵；d是差商矩阵；f是所求自变量对应的因变量；

n=length(x);

d=zeros(n,n);

d(:,1)=y';

f=0;

for j=2:n

for i=j:n

d(i,j)=(d(i,j-1)-d(i-1,j-1))/(x(i)-x(i-j+1)); % 构造差商矩阵 end

end

c=d(n,n);

for k=(n-1):-1:1

c=conv(c,poly(x(k)));

m=length(c);

c(m)=c(m)+d(k,k); % 又插值公式构造各阶次对应系数矩阵

end

c=vpa(c,10); % 显示10位有效数字

for k=1:m

f=f+c(k)*a^(m-k);

end

end

输入矩阵

x=[0 30000 60000 90000 120000]

y=[9.81 9.7487 9.6879 9.6278 9.5682]

a=55000

再运行源函数，可得：

c =

[ -2.057613169e-23, 4.938271605e-18, -3.703703703e-14, -0.000002046111111, 9.81]

d =

9.8100 0 0 0 0

9.7487 -0.0000 0 0 0

9.6879 -0.0000 0.0000 0 0

9.6278 -0.0000 0.0000 0.0000 0

9.5682 -0.0000 0.0000 -0.0000 -0.0000

f =

9.6979851723251030491735014653435

即此时得到的插值函数为：

f(x)=−2.06×10−23x4+4.94×10−18x3−3.70×10−14x2−2.05×10−6x+9.81本题结果为：

f(55000)=9.6980

3.分段低次插值法

当插值点比较多的时候，拉格朗日插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却会和“实际上”的值之间有很大的偏差。这类现象也被称为Runge现象，解决的办法是分段用较低次数的插值多项式。

Runge现象