二元拉格朗日插值Fortran程序设计实验.doc

《程序设计》编程实验 XXX
二元拉格朗日插值

一 实验目的-程序功能
利用FORTRAN编程实现二元拉格朗日插值求解函数在给定点的函数值。
设已知插值节点（xi,yi）(i=1,…,m,j=1,…,n)及对应函数值zij=f(xi,yi) (i=1,…,m,j=1,…,n),用拉格朗日插值法求函数在给定点（x,y）处的对应函数值z。
二 实验内容
1、 了解和学习FORTRAN程序语言，会编写一些小程序；
2、 学习和理解拉格朗日插值的原理及方法，并拓展至二元拉格朗日插值方法；
3、 利用FORTRAN编程实现二元拉格朗日插值法；
4、 举例进行求解，并对结果进行分析。
三 实验原理及方法
1、基本概念
已知函数y=f(x)在若干点的函数值=（i=0,1,,n）一个差值问题就是求一“简单”的函数p(x)：
p()=,i=0,1,,n, (1)
则p(x)为f(x)的插值函数，而f(x)为被插值函数会插值原函数，，，，...,为插值节点，式（1）为插值条件，如果对固定点求f()数值解，我们称为一个插值节点，f()p()称为点的插值，当[min(，，，...,)，max(，，，...,)]时，称为内插，否则称为外插式外推，特别地，当p(x)为不超过n次多项式时称为n阶Lagrange插值。
2、Lagrange插值公式
2.1 线性插值
设已知 ，及=f() ,=f(),为不超过一次多项式且满足=,=,几何上，为过（，），（，）的直线，从而得到
=+（x-）. （2）
为了推广到高阶问题，我们将式（2）变成对称式
=（x）+(x). （3）
其中，
（x）=,(x)=。
均为1次多项式且满足
（）=1且()=0；（）=0且()=1。
（4）
两关系式可统一写成

2.2 n阶Lagrange插值
（5）
设已知，，，...,及=f()(i=0,1,.....,n),为不超过n次多项式且满足（i=0,1,...n）.
易知
其中，均为n次多项式且满足式（4）（i,j=0,1,...,n）,再由（ji）为n次多项式的n个根，知=c.最后，由
c=,i=0,1,...,n.
总之得到：
（6）

（7）

（6）式为n阶Lagrange插值公式，其中，（i=0,1,…,n）称为n阶Lagrange插值的基函数。
3 二元拉格朗日插值方法
对于一元函数y=f(x)，得到n+1个数据点（,）(i=0,1,…,n),可由（6）、（7）式求得n阶Lagrange插值公式，然后求函数在y=f(x)在x点的函数值。
对于二元函数，若知道数据点（i=1,…,m，j=1,…,n）,可利用两次拉格朗日插值计算在点（x,y）的函数值，方法如下：
（1）对每个( i=1,…,m),以（ j=1,…,n）为插值节点，
（ j=1,…,n）为对应函数值，y为插值变量，作一元函数插值得( i=1,…,m);
（2） 以( i=1,…,m)为插值节点，( i=1,…,m)为对应函数值，x为插值变量，作一元函数插值求得（x,y）点的值z。
四 FORTRAN编程
a) 开发环境
使用Compaq Visual Fortran 6.6进行程序设计，编程实现二元拉格朗日插值算法。
b) 使用说明
先编出一元拉格朗日差值算法子程序lagrange，然后编写二元拉格朗日插值算法程序lagrange2,其中两次调用lagrange子程序。
Lagrange(xa,ya,n,x,y)
n 整型变量，输入参数，节点个数
xa n个元素的一维实数型数组，输入参数，存放自变量插值节点xi（i=1,…,n）
ya n个元素的一维实数型数组，输入参数,存放函数值（y1,…,yn）T
x 实型变量，输入参数，插值自变量
y 实型变量，输出参数，所求值
******************************************************
Lagrange2(x1a, x2a,ya,m,n,x1, x2,y)
m 整型变量，输入参数，x自变量节点个数
n 整型变量，输入参数，y自变量节点个数
x1a m个元素的一维实数型数组，输入参数，存放x自变量插值节点xi（i=1,…,m）
x2a n个元素的一维实数型数组，输入参数，存放y自变量插值节点yj（j=1,…,n）
x1 实型变量，输入参数，插值x自变量
x2 实型变量，输入参数，插值y自变量
ya m×n个元素的二维实数型数组，输入参数,存放（xi,yj）（i=1,…,m，j=1,…,n）函数值（y1,…,yn）T
y 实型变量，输出参数，所求插值结果
c) 程序代码
Lagrange子程序
SUBROUTINE lagrange(xa,ya,n,x,y)
integer n,nmax
real x,y,xa(n),ya(n),l(n),dy
parameter(nmax=10)
integer i,j
l(1)=1
do j=2,n
l(1)=l(1)*(x-xa(j))/(xa(1)-xa(j)) !计算l1(x)
end do
do i=2,n-1
l(i)=1
do j=1,i-1
l(i)=l(i)*(x-xa(j))/(xa(i)-xa(j))
end do
do j=i+1,n
l(i)=l(i)*(x-xa(j))/(xa(i)-xa(j)) !计算li(x),1<i<n
end do
end do
l(n)=1
do j=1,n-1
l(n)=l(n)*(x-xa(j))/(xa(n)-xa(j)) !计算ln(x)
end do
y=0
do i=1,n
y=y+l(i)*ya(i) !计算y=求和li(x)*ya(i)
end do
end subroutine lagrange

Lagrange子程序说明：
已知数据点（xa(i),ya(i)）(i=0,1,…,n),求插值多项式,其中
先求，程序中l(n)为一维实型数组，存放插值基函数，l(1)即为;
然后，对于i=2, …,n-1,

最后计算
求和得到 （i=1,2，…,n）
对于每一个自变量x输入参数，可以得到一个y输出参数。
Lagrange2子程序
SUBROUTINE lagrange2(x1a,x2a,ya,m,n,x1,x2,y)
integer n,nmax,m,mmax
real x1,x2,y,x1a(m),x2a(n),ya(m,n)
parameter(nmax=100,mmax=100)
integer i,j
real ym(mmax),yn(nmax)
do i=1,m
do j=1,n
yn(j)=ya(i,j) ！对每一个xi,以（yj,zij）作为插值节点
end do
call lagrange(x2a,yn,n,x2,ym(i)) ！求得（xi,y）的函数值ui
end do
call lagrange(x1a,ym,m,x1,y) !以（xi,ui）插值点求（x,y）函数值
end subroutine lagrange2

Lagrange2子程序说明：
首先对每一个x1=x1a(i)（i=1,2，…,m），也就是x=xi,以（yj,zij）作为插值节点，得到插值多项式，以y为变量，可求得（xi,y）点的函数值ui，程序中调用一次lagrange子程序，以x2a（即为yj，j=1,2，…,n）、yn（即为zij, j=1,2，…,n）输入得到x2=y点(对应点（xi,y）)的值ym(i)(即为ui) （i=1,2，…,m）；
然后以(xi,ui) （i=1,2，…,m）为插值点，得到插值多项式，以x为变量，求得点（x，y）点的函数值z=f(x,y),程序中再次调用lagrange子程序，以x1a（即为xi，i=1,2，…,m）、ym（即为ui, i=1,2，…,m）输入得到x1=x点(对应点（x,y）)的值y,也就是z=f(x，y)使用二元拉格朗日插值法的计算值。
五 举例验证
Lagrange子程序参考了参考书《Visual Basic常用数值算法集》（何光渝，北京航科学出版社，2002）73页～75页，但不相同，参考书中使用了Neville算法，而以上子程序则是使用的拉格朗日插值得基本定义算法。
与参考书75页使用相同的例子进行验证，f(x)=sinx，验证程序如下（见附件la2.f90）：

图1 验证一元拉格朗日算法 f(x)=sinx
program l1
parameter(n=10,pi=3.1415926)
real dy
dimension xa(n),ya(n)
write(*,*)'y=sin(x) 0<x<PI'
write(*,*)'sin function from 0 to PI'
do i=1,n
xa(i)=i*pi/n ! 输入xi
ya(i)=sin(xa(i)) ! 输入yi
end do
write(*,'(T10,A1,T20,A4,T28,A12,T46,A5)')'x','f(x)','interpolated','error'
do i=1,10
x=(-0.05+i/10.0)*pi ! 自变量x
f=sin(x) ! f(x)真实值
call lagrange(xa,ya,n,x,y) !调用子程序lagrange，求解y
dy=y-f !dy为误差，即插值求解值与真实值之差
write(*,'(1x,3F12.6,F15.6)')x,f,y,dy
end do
end program

运行后，得到：

图2 验证f(x)=sinx结果

以上结果与参考书（下图）进行对照

图3 参考书f(x)=sinx验证结果（P76～P77）
通过对比可以看到，误差数量级差不多，说明本子程序是可行的。
同样对f(x)=e^x进行验证，只需将以上程序中的sin更改为exp，变量x进行相应更改（见附件la1.f90）；

图4 f(x)=e^x验证程序

图5 f(x)=e^x验证结果

图6 参考书77页f(x)=e^x验证结果
对比两个验证结果可以看到参考书的插值程序计算的误差比较小（10-11～10-8），而本实验的对f(x)=e^x验证结果误差较大（10-6～10-5，其中误差为0的除外），说明Neville插值方法改善了精度。
下面进行二元拉格朗日插值计算验证，同样实用参考书P104的例子进行验证，函数为f(x,y)=eysinx，0<x<PI，0<y<1。
编写验证程序如下（见附件l2.f90）：
program l2
parameter(n=5,pi=3.1415926)
real dy
dimension x1a(n),x2a(n),ya(n,n)
write(*,*)'f(x)=sin(x)e^y 0<x<PI,0<y<1'
write(*,*)'f(x)=sinx*e^y function from 0 to PI,0 to 1'
do i=1,n
x1a(i)=i*pi/n !输入xi
do j=1,n
x2a(j)=1.0*j/n !输入yj
ya(i,j)=sin(x1a(i))*exp(x2a(j)) !输入ya(i,j),即zij
end do
end do
write(*,'(T10,A1,T22,A,T32,A,T40,A,T58,A)')'x','y','f(x)','interpolated','error'
do i=1,5
x1=(-0.1+i/5.0)*pi !赋予自变量x值
do j=1,5
x2=-0.1+j/5.0 !赋予自变量y值
f=sin(x1)*exp(x2) !f(x,y)真实值
call lagrange2(x1a,x2a,ya,n,n,x1,x2,y) !调用程序lagrange2计算z=f(x,y)
dy=y-f !计算二元拉格朗日插值法的误差
write(*,'(1x,4F12.6,F14.7)')x1,x2,f,y,dy !输出
end do
write(*,*)'*************************************************************'
end do
end program l2

程序中输入数据节点（xi,yj）及函数值zij，调用lagrange2子程序进行求解(x,y)点的函数值z(即为程序中的y)，其中x≠xi，y≠yj，函数在（x，y）处的真实值为f（x，y）（程序中为f），并求解插值法的误差dy=z-f。

图7 f(x)=sin(x)e^y验证程序
运行后得到的验证结果：

图8 二元拉格朗日插值法输出结果

图9 参考书f(x)=sin(x)e^y验证结果
参考书给自变量赋值4×4个点，本实验赋值了5×5个数据点，可看出该实验程序计算的误差值比参考书误差值小，取得了良好的效果。
最后来用几个其他函数进行验证。
1、f(x,y)=(x+y)3
程序见图10（见附件l22.f90），验证结果见图11。

图10 二元拉格朗日插值法计算程序：f(x,y)=(x+y)3

图11 验证结果：f(x,y)=(x+y)3
2、f(x,y)=(x+y5)sin(xy)
程序见图12（见附件l23.f90），验证结果见图13。

图12 二元拉格朗日插值法计算程序：f(x,y)= (x+y5) sin(xy)

图13 验证结果：f(x,y)= (x+y5) sin(xy)
由以上验证结果可以看出用二元拉格朗日插值法计算较简单的函数f(x,y)=(x+y)3值时误差较小（10-6），而用其计算较复杂的函数f(x,y)= (x+y5) sin(xy)时，误差较大（10-3～10-2），这也是符合插值法计算的原理的。
六 实验总结
通过本次的程序语言实验设计，对Fortran程序语言有了一定的认识和了解，能够编写较简单的程序，实现一定的功能；用Fortran编程实现了一元及二元拉格朗日插值法求解函数在某点的值。
学习一种程序语言是充满乐趣和挑战的，在实际工作中将会得到更大的应用，X所XXX研究员给我们介绍了丰富且实用的面向科学计算的软件包，开拓了我们的视野，让我们能认识到还有很多很多好用、实用的知识需要去学习和掌握。

参考文献：
[1] 白云，李学哲，陈国新，贾波等，FORTRAN 95程序设计，清华大学出版社，2011
[2] 何光渝等，Visual Basic常用数值算法集，科学出版社，2002
[3] 颜庆津等，数值分析，北京航空航天大学出版社，1999
[4] 何光渝,高永利等，Visual Fortran常用数值算法集，科学出版社，2002
15