《二元拉格朗日插值Fortran程序设计实验16页》由会员分享，可在线阅读，更多相关《二元拉格朗日插值Fortran程序设计实验16页（16页珍藏版）》请在金锄头文库上搜索。
1、程序设计编程实验 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)在若干点的
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,y
7、,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)-xa
8、(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,x2
9、,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 subrouti
10、ne 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、1415926) 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) !调用子程序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)=sin
