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

本文格式為Word版，下載可任意編輯——二元拉格朗日插值Fortran程序設(shè)計(jì)試驗(yàn)《程序設(shè)計(jì)》編程試驗(yàn)XXX

二元拉格朗日插值

一試驗(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ì)。二試驗(yàn)內(nèi)容

1、了解和學(xué)習(xí)FORTRAN程序語言，會(huì)編寫一些小程序；2、學(xué)習(xí)和理解拉格朗日插值的原理及方法，并拓展至二元拉格

朗日插值方法；

3、利用FORTRAN編程實(shí)現(xiàn)二元拉格朗日插值法；4、舉例進(jìn)行求解，并對(duì)結(jié)果進(jìn)行分析。三試驗(yàn)原理及方法1、基本概念

已知函數(shù)y=f(x)在若干點(diǎn)xi的函數(shù)值yi=f?xi?（i=0,1,???,n）一個(gè)差值問題就是求一“簡(jiǎn)單〞的函數(shù)p(x)：

p(xi)=yi,i=0,1,???,n,(1)

x0，則p(x)為f(x)的插值函數(shù)，而f(x)為被插值函數(shù)會(huì)插值原函數(shù)，

x1，x2，...,xn為插值節(jié)點(diǎn)，式（1）為插值條件，假使對(duì)固定點(diǎn)x求f(x)

??數(shù)值解，我們稱x為一個(gè)插值節(jié)點(diǎn)，f(x)?p(x)稱為x點(diǎn)的插值，當(dāng)

x?[min(x0，x1，x2，...,xn)，max(x0，x1，x2，...,xn)]時(shí)，稱為內(nèi)插，

?????否則稱為外插式外推，特別地，當(dāng)p(x)為不超過n次多項(xiàng)式時(shí)稱為n階Lagrange插值。

0

《程序設(shè)計(jì)》編程試驗(yàn)XXX

2、Lagrange插值公式2.1線性插值L1(1)

設(shè)已知x0，x1及y0=f(x0),y1=f(x1),L1(x)為不超過一次多項(xiàng)式且滿足L1(x0)=y0,L1(x1)=y1,幾何上，L1(x)為過（x0，y0），（x1，y1）的直線，從而得到L1(x)=y0+

y1?y0（x-x0）.（2）x1?x0為了推廣到高階問題，我們將式（2）變成對(duì)稱式

L1(x)=l0（x）y0+l1(x)y1.（3）

其中，

x?x0x?x1l0（x）=x0?x1,l1(x)=x1?x0。

均為1次多項(xiàng)式且滿足

l0（x0）=1且l1(x0)=0；l0（x1）=0且l1(x1)=1。

兩關(guān)系式可統(tǒng)一寫成

?1（i?j）li(xi)???0（i?j）（4）

2.2n階Lagrange插值Ln(x)

設(shè)已知x0，x1，x2，...,xn及yi=f(xi)(i=0,1,,n),Ln(x)為不超過n次多項(xiàng)式且滿足Ln(xi)?yi（i=0,1,...n）.

易知

Ln(x)?l0(x)y0??ln(x)yn??li(x)yi0n（5）

li(x)均為n次多項(xiàng)式且滿足式j(luò)其中，（4）（i,j=0,1,...,n）,再由x（j?i）

為n次多項(xiàng)式li(x)的n個(gè)根，知li(x)=c

?x?xj?0i?inj.最終，由

1

《程序設(shè)計(jì)》編程試驗(yàn)XXX

li(xj)?c?(xi?xj)?1?c=

j?0j?in1?(x?x)ijj?0j?in,i=0,1,...,n.

總之得到：

Ln(x)??li(x)yii?0n

（6）

li(x)??j?0j?inx?xjxi?xj.

（7）

（6）式為n階Lagrange插值公式，其中，階Lagrange插值的基函數(shù)。3二元拉格朗日插值方法

li(x)（i=0,1,?,n）稱為n

對(duì)于一元函數(shù)y=f(x)，得到n+1個(gè)數(shù)據(jù)點(diǎn)（xi,yi）(i=0,1,?,n),可由（6）、（7）式求得n階Lagrange插值公式，然后求函數(shù)在y=f(x)在x點(diǎn)的函數(shù)值。

zi,j?f(xi,yj)z?f(x,y)對(duì)于二元函數(shù)，若知道數(shù)據(jù)點(diǎn)

（i=1,?,m，j=1,?,n）,可利用兩次拉格朗日插值計(jì)算點(diǎn)（x,y）的函數(shù)值，方法如下：

z?f(x,y)在

yjzijxi（1）對(duì)每個(gè)(i=1,?,m),以（j=1,?,n）為插值節(jié)點(diǎn)，

（j=1,?,n）為對(duì)應(yīng)函數(shù)值，y為插值變量，作一元函數(shù)插值得

ui(i=1,?,m);

（2）以

xi(i=1,?,m)為插值節(jié)點(diǎn)，ui(i=1,?,m)為對(duì)應(yīng)函數(shù)值，

x為插值變量，作一元函數(shù)插值求得（x,y）點(diǎn)的值z(mì)。四FORTRAN編程

2

《程序設(shè)計(jì)》編程試驗(yàn)XXX

a)開發(fā)環(huán)境

使用CompaqVisualFortran6.6進(jìn)行程序設(shè)計(jì)，編程實(shí)現(xiàn)二元拉格朗日插值算法。b)使用說明

先編出一元拉格朗日差值算法子程序lagrange，然后編寫二元拉格朗日插值算法程序lagrange2,其中兩次調(diào)用lagrange子程序。

Lagrange(xa,ya,n,x,y)

n整型變量，輸入?yún)?shù)，節(jié)點(diǎn)個(gè)數(shù)

xan個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù)，存放自變量插

值節(jié)點(diǎn)xi（i=1,?,n）

yan個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù),存放函數(shù)值

（y1,?,yn）T

x實(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ù)

x1am個(gè)元素的一維實(shí)數(shù)型數(shù)組，輸入?yún)?shù)，存放x自變

量插值節(jié)點(diǎn)xi（i=1,?,m）

x2an個(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自變量

3

《程序設(shè)計(jì)》編程試驗(yàn)XXX

yam×n個(gè)元素的二維實(shí)數(shù)型數(shù)組，輸入?yún)?shù),存放（xi,yj）

（i=1,?,m，j=1,?,n）函數(shù)值（y1,?,yn）T

y實(shí)型變量，輸出參數(shù)，所求插值結(jié)果c)程序代碼Lagrange子程序

SUBROUTINElagrange(xa,ya,n,x,y)

integern,nmax

realx,y,xa(n),ya(n),l(n),dyparameter(nmax=10)integeri,jl(1)=1doj=2,n

l(1)=l(1)*(x-xa(j))/(xa(1)-xa(j))!計(jì)算l1(x)enddodoi=2,n-1l(i)=1

doj=1,i-1l(i)=l(i)*(x-xa(j))/(xa(i)-xa(j))enddodoj=i+1,n

l(i)=l(i)*(x-xa(j))/(xa(i)-xa(j))!計(jì)算li(x),1《程序設(shè)計(jì)》編程試驗(yàn)XXX

Lagrange子程序說明：

已知數(shù)據(jù)點(diǎn)（xa(i),ya(i)）(i=0,1,?,n),求插值多項(xiàng)式

n,其中n(x?xa(j))先求l1(x)??，程序中l(wèi)(n)為一維實(shí)型數(shù)組，存放

2(xa(1)?xa(j))插值基函數(shù)，l(1)即為l1(x);然后，對(duì)于i=2,?,n-1,

n(x?xa(j))(x?xa(j))l(i)?li(x)?[?]??

j?1(xa(i)?xa(j))j?i?1(xa(i)?xa(j))i?1y??li(x)?ya(i)li(x)??j?0j?ix?xjxi?xj.(x?xa(j))最終計(jì)算l(n)?ln(x)??

(xa(n)?xa(j))j?1求和得到

n?1y??li(x)?ya(i)（i=1,2，?,n）

對(duì)于每一個(gè)自變量x輸入?yún)?shù)，可以得到一個(gè)y輸出參數(shù)。

Lagrange2子程序

SUBROUTINElagrange2(x1a,x2a,ya,m,n,x1,x2,y)

integern,nmax,m,mmax

realx1,x2,y,x1a(m),x2a(n),ya(m,n)parameter(nmax=100,mmax=100)integeri,j

realym(mmax),yn(nmax)doi=1,mdoj=1,n

yn(j)=ya(i,j)！對(duì)每一個(gè)xi,以（yj,zij）作為插值節(jié)點(diǎn)enddo

calllagrange(x2a,yn,n,x2,ym(i))！求得（xi,y）的函數(shù)值uienddo

calllagrange(x1a,ym,m,x1,y)!以（xi,ui）插值點(diǎn)求（x,y）函數(shù)值endsubroutinelagrange2

5

《程序設(shè)計(jì)》編程試驗(yàn)XXX

Lagrange2子程序說明：

首先對(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子程序參考了參考書《VisualBasic常用數(shù)值算法集》（何光渝，北京航科學(xué)出版社，2023）73頁～75頁，但不一致，參考書中使用了Neville算法，而以上子程序則是使用的拉格朗日插值得基本定義算法。

與參考書75頁使用一致的例子進(jìn)行驗(yàn)證，f(x)=sinx