《二元拉格朗日插值Fortran程序设计实验》由会员分享，可在线阅读，更多相关《二元拉格朗日插值Fortran程序设计实验（16页珍藏版）》请在金锄头文库上搜索。
1、程序设计编程实验 XXX0二元拉格朗日插值一 实验目的-程序功能利用 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、基本概念已知函数
2、 y=f(x)在若干点 的函数值 = （i=0,1, ,n）一个ixiyixf差值问题就是求一“简单”的函数 p(x)：p( )= ,i=0,1, ,n, (1)ixiy则 p(x)为 f(x)的插值函数，而 f(x)为被插值函数会插值原函数， ， ，., 为插值节点，式（ 1）为插值条件，如果对固定点0x12n求 f( )数值解，我们称 为一个插值节点，f( ) p( )称为 点的插 xxx值，当 min( ， ， ，., )，max( ， ， ，., )时，称x012n012n为内插，否则称为外插式外推，特别地，当 p(x)为不超过 n 次多项式时称为 n 阶 Lagrange 插值。程序
3、设计编程实验 XXX12、Lagrange 插值公式 2.1 线性插值 )1(L设已知 ， 及 =f( ) , =f( ), 为不超过一次多项式0x0yx1yx)(1L且满足 = , = ,几何上， 为过（ ， ） ， （ ， ）)(1)(1 )(0y1xy的直线，从而得到= + （x- ）. )(1xL0y01x（2）为了推广到高阶问题，我们将式（2）变成对称式= （x） + (x) . )(1L0l0y1l（3）其中，（x）= , (x)= 。0l10xl01x均为 1 次多项式且满足（ ）=1 且 ( )=0； （ ）=0 且 ( )=1。0l1l00l11l两关系式可统一写成i()0l
4、x2.2 n 阶 Lagrange 插值 Ln设已知 ， ， ，., 及 =f( )(i=0,1,.,n), 为不超过 n0x12xiyix)(xLn次多项式且满足 （i=0,1,.n）.iin)(易知 0 0()()()nn niiLxlylxylxyL其中， 均为 n 次多项式且满足式（4） （i,j=0,1,.,n）,再由li（4）j（ ）（ ）（5）程序设计编程实验 XXX2（j i）为 n 次多项式 的 n 个根，知 =c .最后，由x)(xli )(xlinijjx0c= ,i=0,1,.,n.1)()(0ijjji xcl ijj0)(总之得到： 0()()nniiLxly0()
5、.njijijxl（6）式为 n 阶 Lagrange 插值公式，其中， （i=0,1,n）称为)(xlin 阶 Lagrange 插值的基函数。3 二元拉格朗日插值方法对于一元函数 y=f(x)，得到 n+1 个数据点（ , ）(i=0,1,n),ixy可由（6） 、 （7）式求得 n 阶 Lagrange 插值公式，然后求函数在y=f(x)在 x 点的函数值。对于二元函数 ，若知道数据点(,)zfxy（i=1,m，j=1,n）,可利用两次拉格朗日插值计,(,)ijijzfxy算 在点（x,y）的函数值，方法如下：,f（1）对每个 ( i=1,m),以 （ j=1,n）为插值节点，ixjy
6、ijz（ j=1,n）为对应函数值， y 为插值变量，作一元函数插值得 ( iui=1,m);（2） 以 ( i=1,m)为插值节点， ( i=1,m)为对应函数ix iu（6）（7）程序设计编程实验 XXX3值，x 为插值变量，作一元函数插值求得（x,y）点的值 z。四 FORTRAN 编程a) 开发环境使用 Compaq Visual Fortran 6.6 进行程序设计，编程实现二元拉格朗日插值算法。b) 使用说明先编出一元拉格朗日差值算法子程序 lagrange，然后编写二元拉格朗日插值算法程序 lagrange2,其中两次调用 lagrange 子程序。Lagrange(xa,ya,
7、n,x,y)n 整型变量，输入参数，节点个数xa n 个元素的一维实数型数组，输入参数，存放自变量插值节点 xi（i=1,n）ya n 个元素的一维实数型数组，输入参数,存放函数值（y 1,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）
8、程序设计编程实验 XXX4x1 实型变量，输入参数，插值 x 自变量x2 实型变量，输入参数，插值 y 自变量ya mn 个元素的二维实数型数组，输入参数,存放（x i,yj） （i=1,m，j=1,n）函数值（y 1,yn） Ty 实型变量，输出参数，所求插值结果c) 程序代码Lagrange 子程序SUBROUTINE lagrange(xa,ya,n,x,y)integer n,nmaxreal x,y,xa(n),ya(n),l(n),dyparameter(nmax=10)integer i,jl(1)=1do j=2,nl(1)=l(1)*(x-xa(j)/(xa(1)-xa(j)
9、 !计算 l1(x)end dodo i=2,n-1l(i)=1do j=1,i-1l(i)=l(i)*(x-xa(j)/(xa(i)-xa(j) end dodo j=i+1,nl(i)=l(i)*(x-xa(j)/(xa(i)-xa(j) !计算 li(x),1inend doend dol(n)=1do j=1,n-1l(n)=l(n)*(x-xa(j)/(xa(n)-xa(j) !计算 ln(x)end doy=0do i=1,ny=y+l(i)*ya(i) !计算 y=求和 li(x)*ya(i)end doend subroutine lagrange程序设计编程实验 XXX5La
10、grange 子程序说明：已知数据点（xa(i),ya(i)）(i=0,1,n),求插值多项式,其中()(i)iylxya0().njijijxl先求 ，程序中 l(n)为一维实型数组，存12(j)(1nxa放插值基函数，l(1) 即为 ;l然后，对于 i=2, ,n-1, 1 1(j)(j)(i)xaxai nij jillxi i 最后计算1()()()njnnjxll求和得到 （i=1,2，,n）()(i)iylxya对于每一个自变量 x 输入参数，可以得到一个 y 输出参数。Lagrange2 子程序SUBROUTINE lagrange2(x1a,x2a,ya,m,n,x1,x2,y
11、)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,mdo j=1,nyn(j)=ya(i,j) ！对每一个 xi,以（yj,zij ）作为插值节点end docall lagrange(x2a,yn,n,x2,ym(i) ！求得（xi,y ）的函数值 uiend docall lagrange(x1a,ym,m,x1,y) !以（xi,ui）插值点求（x,y）函数值end subroutine la
12、grange2程序设计编程实验 XXX6Lagrange2 子程序说明：首先对每一个 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)
13、,程序中再次调用 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）：程序设计编程实验 XXX7图 1 验证一元拉格朗日算法 f(x)=sinxprogram l1parameter(n=10,pi=3.1415926)real dydimension xa(n),ya(n)write(*,*)y=sin(x) 0xPIwrite(*,*)sin function from 0 to PIdo i=1,nxa(i)=i*pi/n ! 输入 xiya(i)=sin(xa(i) ! 输入 yiend dowrite(*,(T10,A1,T20,A4,T28,A12,T46,A5)x,f(x)
