资源预览内容
第1页 / 共16页
第2页 / 共16页
第3页 / 共16页
第4页 / 共16页
第5页 / 共16页
第6页 / 共16页
第7页 / 共16页
第8页 / 共16页
第9页 / 共16页
第10页 / 共16页
亲,该文档总共16页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述
程序设计编程实验 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)Tx 实型变量,输入参数,插值自变量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 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,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(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,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 subroutine 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(即为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.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
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号