二元拉格朗日插值Fortran程序设计实验(共16页).doc

《二元拉格朗日插值Fortran程序设计实验(共16页).doc》由会员分享，可在线阅读，更多相关《二元拉格朗日插值Fortran程序设计实验(共16页).doc（16页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上二元拉格朗日插值一 实验目的-程序功能利用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)在若干点

2、的函数值=（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

3、）为了推广到高阶问题，我们将式（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 二元拉格朗日插

4、值方法 对于一元函数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

5、 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）Tx 实型变量，输入参数，插值自变量y 实型变量，输出参数，所求值*Lagrange2(x1a, x2a,ya,m,n,x1, x2,y)m 整型变

6、量，输入参数，x自变量节点个数n 整型变量，输入参数，y自变量节点个数x1a m个元素的一维实数型数组，输入参数，存放x自变量插值节点xi（i=1,m）x2a n个元素的一维实数型数组，输入参数，存放y自变量插值节点yj（j=1,n）x1 实型变量，输入参数，插值x自变量x2 实型变量，输入参数，插值y自变量ya mn个元素的二维实数型数组，输入参数,存放（xi,yj）（i=1,m，j=1,n）函数值（y1,yn）Ty 实型变量，输出参数，所求插值结果c) 程序代码Lagrange子程序SUBROUTINE lagrange(xa,ya,n,x,y) integer n,nmaxreal x,

7、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),1in end do end do l(n)=1 do j=1,n-1 l(n)=l(n)*(x-xa(j)/(xa(n)-x

8、a(j) !计算ln(x) end doy=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,x

9、2,y) integer n,nmax,m,mmaxreal x1,x2,y,x1a(m),x2a(n),ya(m,n) parameter(nmax=100,mmax=100)integer i,jreal 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）的函数值uiend do call lagrange(x1a,ym,m,x1,y) !以（xi,ui）插值点求（x,y）函数值end subrout

10、ine lagrange2Lagrange2子程序说明： 首先对每一个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（

11、即为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)=sinxprogram l1 parameter(n=10,pi=3

12、.) real dy dimension xa(n),ya(n) write(*,*)y=sin(x) 0xPI 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) !

13、调用子程序lagrange，求解y dy=y-f !dy为误差，即插值求解值与真实值之差 write(*,(1x,3F12.6,F15.6)x,f,y,dy end doend program运行后，得到：图2 验证f(x)=sinx结果以上结果与参考书（下图）进行对照图3 参考书f(x)=sinx验证结果（P76P77）通过对比可以看到，误差数量级差不多，说明本子程序是可行的。同样对f(x)=ex进行验证，只需将以上程序中的sin更改为exp，变量x进行相应更改（见附件la1.f90）；图4 f(x)=ex验证程序 图5 f(x)=ex验证结果 图6 参考书77页f(x)=ex验证结果对比两

14、个验证结果可以看到参考书的插值程序计算的误差比较小（10-1110-8），而本实验的对f(x)=ex验证结果误差较大（10-610-5，其中误差为0的除外），说明Neville插值方法改善了精度。下面进行二元拉格朗日插值计算验证，同样实用参考书P104的例子进行验证，函数为f(x,y)=eysinx，0xPI，0y1。编写验证程序如下（见附件l2.f90）：program l2 parameter(n=5,pi=3.) real dy dimension x1a(n),x2a(n),ya(n,n) write(*,*)f(x)=sin(x)ey 0xPI,0y1 write(*,*)f(x)=

15、sinx*ey 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(x

16、1)*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 doend program l2程序中输入数据节点（xi,yj）及函数值zij，调用lagrange2子程序进行求解(x,y)点的函数值z(即为程序中的y)，其中xxi，yyj，函数在（x，y）处的真实值为f（x，y）（程序中为f），并求解插值法的误差

17、dy=z-f。 图7 f(x)=sin(x)ey验证程序运行后得到的验证结果：图8 二元拉格朗日插值法输出结果图9 参考书f(x)=sin(x)ey验证结果参考书给自变量赋值44个点，本实验赋值了55个数据点，可看出该实验程序计算的误差值比参考书误差值小，取得了良好的效果。最后来用几个其他函数进行验证。1、f(x,y)=(x+y)3程序见图10（见附件l22.f90），验证结果见图11。 图10 二元拉格朗日插值法计算程序：f(x,y)=(x+y)3图11 验证结果：f(x,y)=(x+y)32、f(x,y)=(x+y5)sin(xy)程序见图12（见附件l23.f90），验证结果见图13。图

18、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-310-2），这也是符合插值法计算的原理的。六 实验总结通过本次的程序语言实验设计，对Fortran程序语言有了一定的认识和了解，能够编写较简单的程序，实现一定的功能；用Fortran编程实现了一元及二元拉格朗日插值法求解函数在某点的值。学习一种程序语言是充满乐趣和挑战的，在实际工作中将会得到更大的应用，X所XXX研究员给我们介绍了丰富且实用的面向科学计算的软件包，开拓了我们的视野，让我们能认识到还有很多很多好用、实用的知识需要去学习和掌握。参考文献：1 白云，李学哲，陈国新，贾波等，FORTRAN 95程序设计，清华大学出版社，20112 何光渝等，Visual Basic常用数值算法集，科学出版社，20023 颜庆津等，数值分析，北京航空航天大学出版社，19994 何光渝,高永利等，Visual Fortran常用数值算法集，科学出版社，2002专心-专注-专业