二元拉格朗日插值Fortran程序設(shè)計(jì)實(shí)驗(yàn)

程序設(shè)計(jì)編程實(shí)驗(yàn) XXX二元拉格朗日插值一 實(shí)驗(yàn)?zāi)康?程序功能利用FORTRAN編程實(shí)現(xiàn)二元拉格朗日插值求解函數(shù)在給定點(diǎn)的函數(shù)值。設(shè)已知插值節(jié)點(diǎn)（xi,yi）(i=1,m,j=1,n)及對(duì)應(yīng)函數(shù)值z(mì)ij=f(xi,yi) (i=1,m,j=1,n),用拉格朗日插值法求函數(shù)在給定點(diǎn)（x,y）處的對(duì)應(yīng)函數(shù)值z(mì)。二 實(shí)驗(yàn)內(nèi)容1、 了解和學(xué)習(xí)FORTRAN程序語(yǔ)言，會(huì)編寫(xiě)一些小程序；2、 學(xué)習(xí)和理解拉格朗日插值的原理及方法，并拓展至二元拉格朗日插值方法；3、 利用FORTRAN編程實(shí)現(xiàn)二元拉格朗日插值法；4、 舉例進(jìn)行求解，并對(duì)結(jié)果進(jìn)行分析。三 實(shí)驗(yàn)原理及方法1、基本概念已知函數(shù)y=f(x)在若干點(diǎn)的函數(shù)值=（i=0,1,n）一個(gè)差值問(wèn)題就是求一“簡(jiǎn)單”的函數(shù)p(x)：p()=,i=0,1,n, (1)則p(x)為f(x)的插值函數(shù)，而f(x)為被插值函數(shù)會(huì)插值原函數(shù)，.,為插值節(jié)點(diǎn)，式（1）為插值條件，如果對(duì)固定點(diǎn)求f()數(shù)值解，我們稱為一個(gè)插值節(jié)點(diǎn)，f()p()稱為點(diǎn)的插值，當(dāng)min(，.,)，max(，.,)時(shí)，稱為內(nèi)插，否則稱為外插式外推，特別地，當(dāng)p(x)為不超過(guò)n次多項(xiàng)式時(shí)稱為n階Lagrange插值。2、Lagrange插值公式 2.1 線性插值設(shè)已知 ，及=f() ,=f(),為不超過(guò)一次多項(xiàng)式且滿足=,=,幾何上，為過(guò)（，），（，）的直線，從而得到 =+（x-）. （2）為了推廣到高階問(wèn)題，我們將式（2）變成對(duì)稱式=（x）+(x). （3）其中，（x）=,(x)=。均為1次多項(xiàng)式且滿足（）=1且()=0；（）=0且()=1。（4）兩關(guān)系式可統(tǒng)一寫(xiě)成 2.2 n階Lagrange插值（5）設(shè)已知，.,及=f()(i=0,1,.,n),為不超過(guò)n次多項(xiàng)式且滿足（i=0,1,.n）.易知 其中，均為n次多項(xiàng)式且滿足式（4）（i,j=0,1,.,n）,再由（ji）為n次多項(xiàng)式的n個(gè)根，知=c.最后，由c=,i=0,1,.,n.總之得到：（6）（7）（6）式為n階Lagrange插值公式，其中，（i=0,1,n）稱為n階Lagrange插值的基函數(shù)。3 二元拉格朗日插值方法 對(duì)于一元函數(shù)y=f(x)，得到n+1個(gè)數(shù)據(jù)點(diǎn)（,）(i=0,1,n),可由（6）、（7）式求得n階Lagrange插值公式，然后求函數(shù)在y=f(x)在x點(diǎn)的函數(shù)值。 對(duì)于二元函數(shù)，若知道數(shù)據(jù)點(diǎn)（i=1,m，j=1,n）,可利用兩次拉格朗日插值計(jì)算在點(diǎn)（x,y）的函數(shù)值，方法如下： （1）對(duì)每個(gè)( i=1,m),以（ j=1,n）為插值節(jié)點(diǎn)，（ j=1,n）為對(duì)應(yīng)函數(shù)值，y為插值變量，作一元函數(shù)插值得( i=1,m);（2） 以( i=1,m)為插值節(jié)點(diǎn)，( i=1,m)為對(duì)應(yīng)函數(shù)值，x為插值變量，作一元函數(shù)插值求得（x,y）點(diǎn)的值z(mì)。四 FORTRAN編程a) 開(kāi)發(fā)環(huán)境使用Compaq Visual Fortran 6.6進(jìn)行程序設(shè)計(jì)，編程實(shí)現(xiàn)二元拉格朗日插值算法。b) 使用說(shuō)明先編出一元拉格朗日差值算法子程序lagrange，然后編寫(xiě)二元拉格朗日插值算法程序lagrange2,其中兩次調(diào)用lagrange子程序。Lagrange(xa,ya,n,x,y)n 整型變量，輸入?yún)?shù)，節(jié)點(diǎn)個(gè)數(shù)xa n個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù)，存放自變量插值節(jié)點(diǎn)xi（i=1,n）ya n個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù),存放函數(shù)值（y1,yn）Tx 實(shí)型變量，輸入?yún)?shù)，插值自變量y 實(shí)型變量，輸出參數(shù)，所求值*Lagrange2(x1a, x2a,ya,m,n,x1, x2,y)m 整型變量，輸入?yún)?shù)，x自變量節(jié)點(diǎn)個(gè)數(shù)n 整型變量，輸入?yún)?shù)，y自變量節(jié)點(diǎn)個(gè)數(shù)x1a m個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù)，存放x自變量插值節(jié)點(diǎn)xi（i=1,m）x2a n個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù)，存放y自變量插值節(jié)點(diǎn)yj（j=1,n）x1 實(shí)型變量，輸入?yún)?shù)，插值x自變量x2 實(shí)型變量，輸入?yún)?shù)，插值y自變量ya mn個(gè)元素的二維實(shí)數(shù)型數(shù)組，輸入?yún)?shù),存放（xi,yj）（i=1,m，j=1,n）函數(shù)值（y1,yn）Ty 實(shí)型變量，輸出參數(shù)，所求插值結(jié)果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) !計(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) !計(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) !計(jì)算ln(x) end doy=0 do i=1,n y=y+l(i)*ya(i) !計(jì)算y=求和li(x)*ya(i) end do end subroutine lagrange Lagrange子程序說(shuō)明： 已知數(shù)據(jù)點(diǎn)（xa(i),ya(i)）(i=0,1,n),求插值多項(xiàng)式,其中 先求，程序中l(wèi)(n)為一維實(shí)型數(shù)組，存放插值基函數(shù)，l(1)即為; 然后，對(duì)于i=2, ,n-1, 最后計(jì)算 求和得到 （i=1,2，,n） 對(duì)于每一個(gè)自變量x輸入?yún)?shù)，可以得到一個(gè)y輸出參數(shù)。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) ！對(duì)每一個(gè)xi,以（yj,zij）作為插值節(jié)點(diǎn) end do call lagrange(x2a,yn,n,x2,ym(i) ！求得（xi,y）的函數(shù)值uiend do call lagrange(x1a,ym,m,x1,y) !以（xi,ui）插值點(diǎn)求（x,y）函數(shù)值end subroutine lagrange2Lagrange2子程序說(shuō)明： 首先對(duì)每一個(gè)x1=x1a(i)（i=1,2，,m），也就是x=xi,以（yj,zij）作為插值節(jié)點(diǎn)，得到插值多項(xiàng)式，以y為變量，可求得（xi,y）點(diǎn)的函數(shù)值ui，程序中調(diào)用一次lagrange子程序，以x2a（即為yj，j=1,2，,n）、yn（即為zij, j=1,2，,n）輸入得到x2=y點(diǎn)(對(duì)應(yīng)點(diǎn)（xi,y）)的值ym(i)(即為ui) （i=1,2，,m）； 然后以(xi,ui) （i=1,2，,m）為插值點(diǎn)，得到插值多項(xiàng)式，以x為變量，求得點(diǎn)（x，y）點(diǎn)的函數(shù)值z(mì)=f(x,y),程序中再次調(diào)用lagrange子程序，以x1a（即為xi，i=1,2，,m）、ym（即為ui, i=1,2，,m）輸入得到x1=x點(diǎn)(對(duì)應(yīng)點(diǎn)（x,y）)的值y,也就是z=f(x，y)使用二元拉格朗日插值法的計(jì)算值。五 舉例驗(yàn)證Lagrange子程序參考了參考書(shū)Visual Basic常用數(shù)值算法集（何光渝，北京航科學(xué)出版社，2002）73頁(yè)75頁(yè)，但不相同，參考書(shū)中使用了Neville算法，而以上子程序則是使用的拉格朗日插值得基本定義算法。與參考書(shū)75頁(yè)使用相同的例子進(jìn)行驗(yàn)證，f(x)=sinx，驗(yàn)證程序如下（見(jiàn)附件la2.f90）： 圖1 驗(yàn)證一元拉格朗日算法 f(x)=sinxprogram l1 parameter(n=10,pi=3.) 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)真實(shí)值 call lagrange(xa,ya,n,x,y) !調(diào)用子程序lagrange，求解y dy=y-f !dy為誤差，即插值求解值與真實(shí)值之差 write(*,(1x,3F12.6,F15.6)x,f,y,dy end doend program運(yùn)行后，得到：圖2 驗(yàn)證f(x)=sinx結(jié)果以上結(jié)果與參考書(shū)（下圖）進(jìn)行對(duì)照?qǐng)D3 參考書(shū)f(x)=sinx驗(yàn)證結(jié)果（P76P77）通過(guò)對(duì)比可以看到，誤差數(shù)量級(jí)差不多，說(shuō)明本子程序是可行的。同樣對(duì)f(x)=ex進(jìn)行驗(yàn)證，只需將以上程序中的sin更改為exp，變量x進(jìn)行相應(yīng)更改（見(jiàn)附件la1.f90）；圖4 f(x)=ex驗(yàn)證程序 圖5 f(x)=ex驗(yàn)證結(jié)果 圖6 參考書(shū)77頁(yè)f(x)=ex驗(yàn)證結(jié)果對(duì)比兩個(gè)驗(yàn)證結(jié)果可以看到參考書(shū)的插值程序計(jì)算的誤差比較?。?0-1110-8），而本實(shí)驗(yàn)的對(duì)f(x)=ex驗(yàn)證結(jié)果誤差較大（10-610-5，其中誤差為0的除外），說(shuō)明Neville插值方法改善了精度。下面進(jìn)行二元拉格朗日插值計(jì)算驗(yàn)證，同樣實(shí)用參考書(shū)P104的例子進(jìn)行驗(yàn)證，函數(shù)為f(x,y)=eysinx，0xPI，0y1。編寫(xiě)驗(yàn)證程序如下（見(jiàn)附件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)=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(x1)*exp(x2) !f(x,y)真實(shí)值 call lagrange2(x1a,x2a,ya,n,n,x1,x2,y) !調(diào)用程序lagrange2計(jì)算z=f(x,y) dy=y-f !計(jì)算二元拉格朗日插值法的誤差 write(*,(1x,4F12.6,F14.7)x1,x2,f,y,dy !輸出 end do write(*,*)* end doend program l2程序中輸入數(shù)據(jù)節(jié)點(diǎn)（xi,yj）及函數(shù)值z(mì)ij，調(diào)用lagrange2子程序進(jìn)行求解(x,y)點(diǎn)的函數(shù)值z(mì)(即為程序中的y)，其中xxi，yyj，函數(shù)在（x，y）處的真實(shí)值為f（x，y）（程序中為f），并求解插值法的誤差dy=z-f。 圖7 f(x)=sin(x)ey驗(yàn)證程序運(yùn)行后得到的驗(yàn)證結(jié)果：圖8 二元拉格朗日插值法輸出結(jié)果圖9 參考書(shū)f(x)=sin(x)ey驗(yàn)證結(jié)果參考書(shū)給自變量賦值44個(gè)點(diǎn)，本實(shí)驗(yàn)賦值了55個(gè)數(shù)據(jù)點(diǎn)，可看出該實(shí)驗(yàn)程序計(jì)算的誤差值比參考書(shū)誤差值小，取得了良好的效果。最后來(lái)用幾個(gè)其他函數(shù)進(jìn)行驗(yàn)證。1、f(x,y)=(x+y)3程序見(jiàn)圖10（見(jiàn)附件l22.f90），驗(yàn)證結(jié)果見(jiàn)圖11。 圖10 二元拉格朗日插值法計(jì)算程序：f(x,y)=(x+y)3圖11 驗(yàn)證結(jié)果：f(x,y)=(x+y)32、f(x,y)=(x+y5)sin(xy)程序見(jiàn)圖12（見(jiàn)附件l23.f90），驗(yàn)證結(jié)果見(jiàn)圖13。圖12 二元拉格朗日插值法計(jì)算程序：f(x,y)= (x+y5) sin(xy)圖13 驗(yàn)證結(jié)果：f(x,y)= (x+y5) sin(xy)由以上驗(yàn)證結(jié)果可以看出用二元拉格朗日插值法計(jì)算較簡(jiǎn)單的函數(shù)f(x,y)=(x+y)3值時(shí)誤差較?。?0-6），而用其計(jì)算較復(fù)雜的函數(shù)f(x,y)= (x+y5) sin(xy)時(shí)，誤差較大（10-310-2），這也是符合插值法計(jì)算的原理的。六 實(shí)驗(yàn)總結(jié)通過(guò)本次的程序語(yǔ)言實(shí)驗(yàn)設(shè)計(jì)，對(duì)Fortran程序語(yǔ)言有了一定的認(rèn)識(shí)和了解，能夠編寫(xiě)較簡(jiǎn)單的程序，實(shí)現(xiàn)一定的功能；用Fortran編程實(shí)現(xiàn)了一元及