二元拉格朗日插值Fortran程序設(shè)計實驗

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