北師大數(shù)值分析作業(yè)(一)、(二)、(三)==

1、北師大數(shù)值分析作業(yè)（一）、（二）、（三）（一）1、 設(shè)下列各數(shù)均為經(jīng)過四舍五入后得到的近似值，試求各數(shù)的絕對誤差限和相對誤差限。 2、 已知是經(jīng)過四舍五入后得到的近似值，問有幾位有效數(shù)字？3、 計算球的體積，為使其相對誤差限為1%，測量半徑R時，相對誤差最大為多少？1、 分別用Gauss消去法、列主元素法和全主元素法解下列方程組，計算過程保留3位小數(shù)。   2、 用三角分解法求解題1中的方程組。 3、 用緊湊格式解下列方程組，并寫出L,U矩陣。 4、 若求證：(1) (2)    5、 用三角分

2、解法求下列矩陣的逆矩陣。6、 設(shè)有方程組Ax=b，其中 （1） （1）    求出A能進(jìn)行Cholesky分解，即A=LLT（其中L為下三角矩陣）的a取值范圍。（2） （2）    取a=1，對矩陣A進(jìn)行Cholesky分解，并用平方根法求解上述方程組，計算過程保留2位小數(shù)。7、 用追趕法解下列方程組8、 已知 求及，并說明方程組Ax=b是否病態(tài)。9、 已知方程組的解為 (1)計算系數(shù)矩陣的條件數(shù)。(2)取，分別計算殘量。 10、 求解超定方程組的最小二乘解。1、已知函數(shù)表為 -1010512（1） （1

3、）  利用線性插值計算的近似值并估計誤差。（2） （2）  利用二次插值計算的近似值并估計誤差。 2、已知函數(shù)表為 0527072708070927001075001219001188001426 用二次插值計算的近似值。3、已知函數(shù)表為 1346-75814試求其3階Lagrange插值多項式，并以此計算f(2)的近似值。 4、設(shè)為n+1個互異節(jié)點，為這組節(jié)點上的n次Lagrange插值基函數(shù)，試證：   5、已知函數(shù)表為 1615163417021828241450246259265271303035試求其3階

4、Newton插值多項式，并以此計算f(1.682)的近似值。  6、已知函數(shù)表為 123435915分別用Newton向前、向后插值公式計算f(1.5)，f(3.7)的近似值。7、設(shè)，求差商 8、設(shè)是一個n次多項式，試證： 9、設(shè)節(jié)點，求的3次Hermite插值多項式及f(1.03)的近似值，并估計誤差。（二）計算實習(xí)說明書目的：訓(xùn)練運用計算機(jī)進(jìn)行科學(xué)與工程計算的能力。要求：1獨立進(jìn)行算法設(shè)計、程序設(shè)計和上機(jī)運算，并得出正確的結(jié)果。2編制程序時全部采用雙精度，要求按題目的要求設(shè)計輸出，并執(zhí)行打印。3只能根據(jù)題目給出的信息并且只允許一次計算得出全部結(jié)果。題

5、目：第二題 關(guān)于，的下列方程組以及關(guān)于，的下列二維數(shù)表ztu00.40.81.21.62.00-0.5-0.340.140.942.063.50.2-0.42-0.5-0.260.31.182.380.4-0.18-0.5-0.5-0.180.461.420.60.22-0.34-0.58-0.5-0.10.620.80.78-0.020.5-0.66-0.5-0.021.01.50.46-0.26-0.66-0.74-0.5確定了一個二元函數(shù)。1.試用數(shù)值方法求出在區(qū)域上的一個近似表達(dá)式：要求一最小的值達(dá)到以下精度：其中。2.計算的值，以觀察逼近的效果，其中。說明：1用迭代方法求解非線性方程

6、組時，要求近似解向量滿足：2作二元插值時，要使用分片二次代數(shù)插值。3要求由程序自動確定最小的值。4打印以下內(nèi)容： （1）算法的設(shè)計方案。（2）全部源程序（要求注明主程序和每個子程序的功能）。（3）數(shù)表：。（4）選擇過程的值。（5）達(dá)到精度要求時的值以及中的系數(shù)（6）數(shù)表：。5采用f型輸出的準(zhǔn)確值，其余實型采用e型輸出并至少顯示12位有效數(shù)字。一、程序算法的設(shè)計根據(jù)題目要求，并結(jié)合題目的特點，算法設(shè)計方案如下表所示：程序算法的設(shè)計方案（流程圖）利用Doolittle法求解線性方程組。求矩陣F(x(k)以及向量F(x(k)。初始化x，y。及存放的數(shù)組A。結(jié)束輸出x*,y*,z*,p*的數(shù)表并比較逼

7、近的效果。利用的值，進(jìn)行曲面擬合，輸出過程中k,sigma值，滿足精度時，輸出k,sigma值以及擬合系數(shù)重新取*,*利用牛頓法和分片二次插值求出*的值。利用已經(jīng)求出的擬合系數(shù)計算近似表達(dá)式的值*。利用和的值進(jìn)行分片二次插值，得到的近似值，并輸出的數(shù)表。利用牛頓法求解非線性方程組，可以解得：的值開始對應(yīng)的子程序插值定位子程序確定插值的位置。利用分片二次插值求和對應(yīng)的。曲面擬合子程序來確定逼近多項式。求擬合系數(shù)子程序，確定各擬合系數(shù)。矩陣的乘法、轉(zhuǎn)置及求擬子程序。說明：1、此流程圖只是大概概括程序的流程，有些細(xì)節(jié)未能全寫在上面。2、對于求解擬合系數(shù)子程序中采用了矩陣的求逆，主要考慮到求逆的矩陣為

8、方陣，求逆運算和采用高斯消元法原理一樣，求逆顯得更簡潔。3、由于Fortran中數(shù)組的下標(biāo)默認(rèn)都是從1開始的，所以在有些地方造成不便，所以有些地方有加1和減1現(xiàn)象。4、本人使用的編譯器為Fortran Powerstation 4.0，在雙精度數(shù)據(jù)方面存在一些問題，主要是直接輸入一個實數(shù)，其存放的數(shù)據(jù)與原數(shù)存在左右的誤差，所以實數(shù)的輸入中都在后面加了D+n，n為次數(shù)。二、全部源程序本人編程所用的軟件是Fortran，所編寫的程序如下：!主程序! 功能：控制整個程序的運行 ! program mainreal(8),parameter:JD=1.0D-12integer i,j,k,H,Lrea

9、l(8) Temp,A(4),F(4),Deta(4),Jac(4,4),Xn(11),Yn(21)real(8),dimension(11,21):t,u,Zreal(8),external:Fanshuopen(98,file="計算結(jié)果.txt",status="unknown")call InitialXY(Xn,Yn) !初始化X,Y向量do i=1,11 do j=1,21do k=1,4 A(k)=1.0D+0 !初始化迭代向量 deta(k)=0.0D+0end doTemp=1do while(Temp.GT.JD) !循環(huán)結(jié)構(gòu)體 ca

10、ll CalF(A,Xn(i),Yn(j),F) !求F向量 do k=1,4 F(k)=-F(k) end do call Jacobi(A,Jac) !求系數(shù)矩陣 call Doolittle(Jac,F,deta) !求deta Temp=abs(Fanshu(deta)/Fanshu(A)!求精度 do k=1,4 A(k)=A(k)+deta(k) !求新的迭代向量 end doend do t(i,j)=A(1) !把滿足精度的向量 u(i,j)=A(2) !分別存在對應(yīng)向量中 end doend dowrite(98,*)" 計算結(jié)果如下 "write(98,

11、*)"="do i=1,11 do j=1,21 call locate(t(i,j),u(i,j),H,L) !插值的定位call CalZ(t(i,j),u(i,j),H,L,Z(i,j) !插值求函數(shù) write(98,80)"X(",i,")=",Xn(i),",Y(",j,")=",Yn(j) write(98,100)"f(",i,",",j,")=",Z(i,j) write(98,*)"=" end

12、doend docall Nihe(Xn,Yn,Z) !曲面擬合80 format(A3,I2,A2,F9.7,A3,I2,A2,F9.7)100 format(1X,A2,I2,A1,I2,A2,E18.12)end!主程序結(jié)束!以下為子程序!初始化X，Y向量!subroutine InitialXY(Xn,Yn)integer ireal(8) Xn(11),Yn(21)do i=1,11 Xn(i)=(.8D-1)*(i-1)end dodo i=1,21 Yn(i)=(.5D+0)+(.5D-1)*(i-1)end doend subroutine!非線性方程組!subroutine

13、CalF(A,x,y,F)real(8) x,yreal(8),dimension(4):A,FF(1)=(.5D+0)*cos(A(1)+A(2)+A(3)+A(4)-x-(.267D+1)F(2)=A(1)+(.5D+0)*sin(A(2)+A(3)+A(4)-y-(.107D+1)F(3)=(.5D+0)*A(1)+A(2)+cos(A(3)+A(4)-x-(.374D+1)F(4)=A(1)+(.5D+0)*A(2)+A(3)+sin(A(4)-y-(.79D+0)end subroutine!一階導(dǎo)數(shù)矩陣矩陣!subroutine Jacobi(A,Jac) real(8) A(4)

14、,Jac(4,4) Jac(1,1)=-(.5D+0)*sin(A(1) Jac(1,2)=1D+0Jac(1,3)=1D+0Jac(1,4)=1D+0 Jac(2,1)=1D+0Jac(2,2)=(.5D+0)*cos(A(2) Jac(2,3)=1D+0Jac(2,4)=1D+0 Jac(3,1)=0.5D+0Jac(3,2)=1D+0Jac(3,3)=-sin(A(3) Jac(3,4)=1D+0 Jac(4,1)=1D+0Jac(4,2)=0.5D+0Jac(4,3)=1D+0Jac(4,4)=cos(A(4) end subroutine!求無窮范數(shù)!Function Fanshu(

15、S)integer ireal(8) Fanshu,S(4)Fanshu=0.0do i=1,4 if(abs(S(i).GT.Fanshu)then Fanshu=S(i) end ifend doend function!杜利特法解線性方程組!Subroutine Doolittle(A,B,X) integer i,k,j,treal(8) Y,A(4,4),B(4),X(4) do i=2,4 A(i,1)=A(i,1)/A(1,1) end do do k=2,4 do j=k,4 do t=1,k-1 Y=Y+A(k,t)*A(t,j)end doA(k,j)=A(k,j)-Y Y

16、=0.0D+0 end do if(A(k,k).EQ.0) THENwrite(*,*)'錯誤!' end if do i=k+1,4do t=1,k-1 Y=Y+A(i,t)*A(t,k)end doA(i,k)=(A(i,k)-Y)/A(k,k)Y=0.0D+0 end do end do do i=2,4 do t=1,i-1Y=Y+A(i,t)*B(t) end do B(i)=B(i)-Y Y=0.0D+0end do X(4)=B(4)/A(4,4)do i=1,3 j=4-i do t=j+1,4Y=Y+A(j,t)*X(t) end do X(j)=(B(j)

17、-Y)/A(j,j) Y=0.0D+0end do end subroutine!插值的定位! subroutine locate(t,u,H,L)integer i,L,Hreal(8) t,udo i=1,6 if(t.LE.(.1D+0)then H=1 else if(t.GT.(.9D+0)then H=5 else if(t.GT.(.2D+0)*i-(.1D+0).AND.t.LE.(.2D+0)*i+(.1D+0)then H=i+1 end if if(u.LE.(.2D+0)then L=1 else if(u.GT.(1.8D+0)then L=5 else if(u.G

18、T.(.4D+0)*i-(.2D+0).AND.u.LE.(.4D+0)*i+(.2D+0)then L=i+1 end ifend doend subroutine!利用分片二次插值求函數(shù)值! subroutine CalZ(t,u,H,L,Z)integer i,j,k,H,Lreal(8) t,u,Lk,Lr,Zreal(8),dimension(6):t_old,u_oldreal(8),dimension(6,6):z_olddo i=1,6 t_old(i)=(.2D+0)*(i-1) u_old(i)=(.4D+0)*(i-1)end do Data z_old /-.5D+0,

19、-.42D+0,-.18D+0,.22D+0,.78D+0,.15D+1,& -.34D+0,-.5D+0,-.5D+0,-.34D+0,-.2D-1,.46D+0,& .14D+0,-.26D+0,-.5D+0,-.58D+0,-.5D+0,-.26D+0,& .94D+0,.3D+0,-.18D+0,-.5D+0,-.66D+0,-.66D+0,& .206D+1,.118D+1,.46D+0,-.1D+0,-.5D+0,-.74D+0,& .35D+1,.238D+1,.142D+1,.62D+0,-.2D-1,-.5D+0/Z=0.0D+0do

20、i=H-1,H+1 do j=L-1,L+1 Lk=1.0D+0 do k=H-1,H+1 if(k.NE.i)then Lk=Lk*(t-t_old(k)/(t_old(i)-t_old(k) end ifend doLr=1.0D+0do k=L-1,L+1 if(k.NE.j)then Lr=Lr*(u-u_old(k)/(u_old(j)-u_old(k) end ifend doZ=Z+Lk*Lr*z_old(i,j) end doend doend subroutine!矩陣的乘法運算!subroutine MultJZ(A,B,C,m,n,nn)integer m,n,nn,i,

21、j,kreal(8):A(m,nn),B(m,n),C(n,nn)do i=1,m do j=1,nn A(i,j)=0.0D+0do k=1,n A(i,j)=A(i,j)+B(i,k)*C(k,j)end do end doend doend subroutine!矩陣的轉(zhuǎn)置運算!subroutine Zhuanzhi(AT,A,N,M)integer i,j,M,Nreal(8):A(N,M),AT(M,N)do i=1,M do j=1,N AT(i,j)=A(j,i) end doend doend subroutine!矩陣的求逆運算!subroutine Qiuni(AN,A,N

22、)integer i,j,Nreal(8):A(N,N),AN(N,N),B(N,N)do i=1,N do j=1,N if(i.EQ.j)then AN(i,j)=1.0D+0else AN(i,j)=0.0D+0end if end doend doB=Acall Upper(B,AN,N)call Lower(B,AN,N)do i=1,N do j=1,N AN(i,j)=AN(i,j)/B(i,i) end doend doend subroutine!矩陣上三角化程序!subroutine Upper(M,S,N)integer Nreal(8) M(N,N),S(N,N),Ei

23、nteger i,jdo i=1,N do j=i+1,N E=M(j,i)/M(i,i)M(j,i:N)=M(j,i:N)-M(i,i:N)*ES(j,:)=S(j,:)-S(i,:)*E end doend doend subroutine!上三角矩陣對角化程序!subroutine Lower(M,S,N)integer Nreal(8) M(N,N),S(N,N),Einteger i,j do i=N,2,-1 do j=i-1,1,-1 E=M(j,i)/M(i,i)M(j,1:N)=M(j,1:N)-M(i,1:N)*ES(j,:)=S(j,:)-S(i,:)*E end doe

24、nd do end subroutine!求擬合系數(shù)矩陣! subroutine CalCRS(Xn,Yn,k,F,CRS)integer i,j,kreal(8):Xn(11),Yn(21),F(11,21),B(11,k),BT(k,11),Temp1(k,11),&G(21,k),GT(k,21),Temp2(21,k),Temp3(k,21)real(8),dimension(k,k):CRS,BTB,GTG,BTBN,GTGNdo i=1,11 do j=1,k B(i,j)=Xn(i)*(j-1) end doend dodo i=1,21 do j=1,k G(i,j)=

25、Yn(i)*(j-1) end doend docall Zhuanzhi(BT,B,11,k)call Zhuanzhi(GT,G,21,k)call MultJZ(BTB,BT,B,k,11,k) !call Qiuni(BTBN,BTB,k)call MultJZ(GTG,GT,G,k,21,k)! call Qiuni(GTGN,GTG,k) call multJZ(Temp1,BTBN,BT,k,k,11) !call multJZ(Temp2,G,GTGN,21,k,k) !call multJZ(Temp3,Temp1,F,k,11,21) !call multJZ(CRS,Te

26、mp3,Temp2,k,21,k)end subroutine!最小二乘曲面擬合! subroutine Nihe(Xn,Yn,F)integer:i,j,k,r,s,H,Lreal(8):Xn(11),Yn(21),sigma,P(11,21),F(11,21),A(4),deta(4),Tempreal(8):FF(4),Jac(4,4),t(11,21),u(11,21),Z(11,21)real(8),allocatable:C(:,:),CRS(:,:)real(8),external:Fanshu k=0sigma=1write(98,*)"=過程中k和sigma如下=&

27、quot;do while(sigma.GT.1.0E-7.AND.k.LE.10)allocate(C(k+1,k+1)call CalCRS(Xn,Yn,k+1,F,C)write(98,'(1X,A2,I1)')'k=',kdo i=1,11 do j=1,21 P(i,j)=0.0D+0do r=1,k+1 do s=1,k+1 P(i,j)=P(i,j)+C(r,s)*(Xn(i)*(r-1)*(Yn(j)*(s-1) end doend do end doend dosigma=0.0D+0do i=1,11 do j=1,21 sigma=sigm

28、a+(F(i,j)-P(i,j)*2 end doend dowrite(98,180)'signa=',sigmawrite(98,*)"="deallocate(C)k=k+1end doallocate(CRS(k,k)call CalCRS(Xn,Yn,k,F,CRS)write(98,*)"=最終的k和sigma如下="write(98,'(1X,A2,I1)')'k=',k-1write(98,180)'sigma=',sigmawrite(98,*)'=最終的擬合系數(shù)C

29、(r,s)如下='do i=1,k do j=1,k write(98,200)'C(',i-1,',',j-1,')=',CRS(i,j) end doend docall InitialXY2(Xn,Yn)do i=1,11 do j=1,21do k=1,4 A(k)=1.0D+0 deta(k)=0.0D+0end doTemp=1do while(Temp.GT.1.0D-12) call CalF(A,Xn(i),Yn(j),FF) do k=1,4 FF(k)=-FF(k) end do call Jacobi(A,Jac)

30、 call Doolittle(Jac,FF,deta) Temp=abs(Fanshu(deta)/Fanshu(A) do k=1,4 A(k)=A(k)+deta(k) end doend do t(i,j)=A(1) u(i,j)=A(2) end doend dowrite(98,*)"=逼近結(jié)果如下="do i=1,11 do j=1,21 call locate(t(i,j),u(i,j),H,L)call CalZ(t(i,j),u(i,j),H,L,Z(i,j) end doend do do i=1,8 do j=1,5 P(i,j)=0.0D+0do

31、r=1,6 do s=1,6 P(i,j)=P(i,j)+CRS(r,s)*(Xn(i)*(r-1)*(Yn(j)*(s-1) end doend do write(98,160)"X(",i,")=",Xn(i),",Y(",j,")=",Yn(j) write(98,200)"f(",i,",",j,")=",Z(i,j)write(98,200)'P(',i,',',j,')=',P(i,j) writ

32、e(98,*)"=" end doend do160 format(1X,A3,I1,A2,F9.7,A4,I1,A2,F9.7)180 format(1X,A6,E18.12)200 format(1X,A3,I1,A1,I1,A2,E18.12)end subroutine!重新初始化X,Y向量!subroutine InitialXY2(Xn,Yn)integer ireal(8) Xn(11),Yn(21)do i=1,11 Xn(i)=(.1D+0)*iend dodo i=1,21 Yn(i)=(.5D+0)+(.2D+0)*iend doend subrout

33、ine!三、計算結(jié)果以下為“計算結(jié)果.txt”文件中的內(nèi)容，即為計算結(jié)果： 計算結(jié)果如下 = X( 1)= .0000000,Y( 1)= .5000000 f( 1, 1)= .446504018480E+00 = X( 1)= .0000000,Y( 2)= .5500000 f( 1, 2)= .324683262927E+00 = X( 1)= .0000000,Y( 3)= .6000000 f( 1, 3)= .210159686683E+00 = X( 1)= .0000000,Y( 4)= .6500000 f( 1, 4)= .103043608316E+00 = X( 1)

34、= .0000000,Y( 5)= .7000000 f( 1, 5)= .340189556266E-02 = X( 1)= .0000000,Y( 6)= .7500000 f( 1, 6)=-.887358136380E-01 = X( 1)= .0000000,Y( 7)= .8000000 f( 1, 7)=-.173371632750E+00 = X( 1)= .0000000,Y( 8)= .8500000 f( 1, 8)=-.250534611467E+00 = X( 1)= .0000000,Y( 9)= .9000000 f( 1, 9)=-.320276506388E+

35、00 = X( 1)= .0000000,Y(10)= .9500000 f( 1,10)=-.382668066110E+00 = X( 1)= .0000000,Y(11)=1.0000000 f( 1,11)=-.437795766738E+00 = X( 1)= .0000000,Y(12)=1.0500000 f( 1,12)=-.485758941444E+00 = X( 1)= .0000000,Y(13)=1.1000000 f( 1,13)=-.526667254884E+00 = X( 1)= .0000000,Y(14)=1.1500000 f( 1,14)=-.5606

36、38479797E+00 = X( 1)= .0000000,Y(15)=1.2000000 f( 1,15)=-.587796538768E+00 = X( 1)= .0000000,Y(16)=1.2500000 f( 1,16)=-.608269779090E+00 = X( 1)= .0000000,Y(17)=1.3000000 f( 1,17)=-.622189452876E+00 = X( 1)= .0000000,Y(18)=1.3500000 f( 1,18)=-.629688378186E+00 = X( 1)= .0000000,Y(19)=1.4000000 f( 1,

37、19)=-.630899760003E+00 = X( 1)= .0000000,Y(20)=1.4500000 f( 1,20)=-.625956152545E+00 = X( 1)= .0000000,Y(21)=1.5000000 f( 1,21)=-.614988546609E+00 = X( 2)= .0800000,Y( 1)= .5000000 f( 2, 1)= .638015226510E+00 = X( 2)= .0800000,Y( 2)= .5500000 f( 2, 2)= .506611755146E+00 = X( 2)= .0800000,Y( 3)= .600

38、0000 f( 2, 3)= .382176369277E+00 = X( 2)= .0800000,Y( 4)= .6500000 f( 2, 4)= .264863491154E+00 = X( 2)= .0800000,Y( 5)= .7000000 f( 2, 5)= .154780200285E+00 = X( 2)= .0800000,Y( 6)= .7500000 f( 2, 6)= .519926834909E-01 = X( 2)= .0800000,Y( 7)= .8000000 f( 2, 7)=-.434680402049E-01 = X( 2)= .0800000,Y

39、( 8)= .8500000 f( 2, 8)=-.131601056789E+00 = X( 2)= .0800000,Y( 9)= .9000000 f( 2, 9)=-.212431088309E+00 = X( 2)= .0800000,Y(10)= .9500000 f( 2,10)=-.286004551058E+00 = X( 2)= .0800000,Y(11)=1.0000000 f( 2,11)=-.352386078979E+00 = X( 2)= .0800000,Y(12)=1.0500000 f( 2,12)=-.411655456522E+00 = X( 2)=

40、.0800000,Y(13)=1.1000000 f( 2,13)=-.463904911519E+00 = X( 2)= .0800000,Y(14)=1.1500000 f( 2,14)=-.509236724701E+00 = X( 2)= .0800000,Y(15)=1.2000000 f( 2,15)=-.547761117962E+00 = X( 2)= .0800000,Y(16)=1.2500000 f( 2,16)=-.579594388339E+00 = X( 2)= .0800000,Y(17)=1.3000000 f( 2,17)=-.604857258890E+00

41、 = X( 2)= .0800000,Y(18)=1.3500000 f( 2,18)=-.623673421332E+00 = X( 2)= .0800000,Y(19)=1.4000000 f( 2,19)=-.636168248413E+00 = X( 2)= .0800000,Y(20)=1.4500000 f( 2,20)=-.642467656690E+00 = X( 2)= .0800000,Y(21)=1.5000000 f( 2,21)=-.642697102700E+00 = X( 3)= .1600000,Y( 1)= .5000000 f( 3, 1)= .840081395765E+00 = X( 3)= .1600000,Y( 2)= .5500000 f( 3, 2)= .699764165673E+00 = X( 3)= .1600000,Y( 3)= .6000000 f( 3, 3)= .56