北航数值分析大作业3

《数值分析B》

第三次数值分析大作业

院系：04 能源与动力工程学院 姓名： 王 开 逍 学号： SY1104207

一．算法设计方案：

1、使用牛顿迭代法，对原题中给出的X(I)=0.08*I，Y(J)=0.5+0.05*J的11*21组X(I),Y(J)分别求得原题非线性方程组的一组解，于是得到一对和X(I),Y(J)对应的T(I),U(J).

2、对于已经求出来的U(I)，T(J)，使用分片二次代数插值法对原题中关于Z,T,U的数表进行插值，得到Z(I,J).于是产生了Z=F(X,Y)的11*21个数值解.

3、从K=1开始逐渐增大K的值，并使用最小二乘曲面拟合法对Z=F(X,Y)进行拟合，得到每次的K和精度TAO，当TAO

4、计算F(X*,Y*),P(X*,Y*)的值，以观察P(X,Y)逼近F(X,Y)的效果，其中X*(I)=0.1*I,Y*(J)=0.5+0.2*J.

二．源程序：

该计算程序使用FORTRAN语言编写

module linear3 !创建模块LINEAR3,用于被主程序调用 use imsl

implicit none

contains

real(kind=8) function max12(a,b,c,d) !求向量的无穷范数的子程序 real(kind=8) a,b,c,d real(kind=8) a0,b0,c0,d0 a0=abs(a) b0=abs(b) c0=abs(c) d0=abs(d) if(a0>b0)then b0=a0 end if if(b0>c0)then c0=b0 end if if(c0>d0)then d0=c0 end if

max12=d0

end function

subroutine qiugeng(x,y,t,u) !求题中所给非线性方程组的解，并且使用牛顿迭代法 real(kind=8) x,y,t,u

real(kind=8) w,v,a(4,4),det(4),b(4) integer i,j t=0.5 u=0.5 w=1 v=1

do while(1>0) a(1,1)=0-0.5*sin(t) a(2,2)=0.5*cos(u) a(3,3)=0-sin(v) a(4,4)=cos(w) do i=1,4 do j=1,4 if(i/=j)then a(i,j)=1 end if end do end do a(3,1)=0.5 a(4,2)=0.5

b(1)=2.67-0.5*cos(t)-u-v-w+x b(2)=1.07-t-0.5*sin(u)-v-w+y b(3)=3.74-0.5*t-u-cos(v)-w+x b(4)=0.79-t-0.5*u-v-sin(w)+y det=a.ix.b

if(max12(det(1),det(2),det(3),det(4))/max12(t,u,v,w)<=0.1e-11) exit t=t+det(1) u=u+det(2) v=v+det(3) w=w+det(4) end do

end subroutine qiugeng

real(kind=8) function p22(x,y,x0,y0,h,tao,z) !使用分片二次插值法，对Z=F(X,Y)进行插值 real(kind=8) x,y,x0(:),y0(:),h,tao,z(:,:) real(kind=8) lx(3),ly(3) integer i,j,n,m,k1,k2 m=size(y0,1) n=size(x0,1)

do i=3,n-2

if(x>(x0(i)-0.5*h) .and. x<=(x0(i)+0.5*h)) then k1=i end if end do

if(x<=(x0(2)+0.5*h)) then k1=2

end if

if(x>(x0(n-1)-0.5*h)) then k1=n-1

end if

do i=3,m-2

if(y>(y0(i)-0.5*tao) .and. y<=(y0(i)+0.5*tao)) then k2=i end if end do

if(y<=(y0(2)+0.5*tao)) then k2=2 end if

if(y>(y0(m-1)-0.5*tao)) then k2=m-1 end if

lx(1)=(x-x0(k1))*(x-x0(k1+1))/((x0(k1-1)-x0(k1))*(x0(k1-1)-x0(k1+1))) lx(2)=(x-x0(k1-1))*(x-x0(k1+1))/((x0(k1)-x0(k1-1))*(x0(k1)-x0(k1+1))) lx(3)=(x-x0(k1))*(x-x0(k1-1))/((x0(k1+1)-x0(k1))*(x0(k1+1)-x0(k1-1))) ly(1)=(y-y0(k2))*(y-y0(k2+1))/((y0(k2-1)-y0(k2))*(y0(k2-1)-y0(k2+1))) ly(2)=(y-y0(k2-1))*(y-y0(k2+1))/((y0(k2)-y0(k2-1))*(y0(k2)-y0(k2+1))) ly(3)=(y-y0(k2-1))*(y-y0(k2))/((y0(k2+1)-y0(k2))*(y0(k2+1)-y0(k2-1))) p22=0.0 do i=1,3 do j=1,3

p22=p22+lx(i)*ly(j)*z(k2-2+j,k1-2+i) end do end do

end function

real(kind=8) function pxy(x,y,c) !关于F(X,Y)的X,Y近似多项式表达式的值 real(kind=8) x,y,c(:,:) integer i,j,m,n real(kind=8) t,p pxy=0.0 m=size(c,1) n=size(c,2) do i=1,m do j=1,n t=x**(i-1) p=y**(j-1) pxy=pxy+c(i,j)*t*p end do end do end function

end module

program main !主程序开始的地方

use linear3 !调用上面的LINEAR3模块

use imsl !调用IMSL库函数，用于曲面拟合的矩阵计算 implicit none

real(kind=8),allocatable:: u(:,:),t(:,:),x(:),y(:),u0(:),t0(:),f(:,:) real(kind=8),allocatable:: b(:,:),g(:,:),bt(:,:),gt(:,:),btb(:,:),gtg(:,:),c(:,:) real(kind=8),allocatable:: btuu(:,:),uu(:,:),btbni(:,:),gtgni(:,:),bniu(:,:),ggni(:,:) real(kind=8) z(6,6),xn(8),yn(5) real(kind=8) h,tao,jingdu integer i,j,k

allocate(u(11,21),t(11,21),x(11),y(21),u0(6),t0(6),f(11,21))

data z /-0.5,-0.42,-0.18,0.22,0.78,1.5,-0.34,-0.5,-0.5,-0.34,-0.02,0.46,& &0.14,-0.26,-0.5,-0.58,-0.5,-0.26,0.94,0.3,-0.18,-0.5,-0.66,-0.66,&

&2.06,1.18,0.46,-0.1,-0.5,-0.74,3.5,2.38,1.42,0.62,-0.02,-0.5/

open(unit=30,file='output/thicknessandf.txt') !将程序运行的结果存储在一个TXT的文档里面

h=0.4 tao=0.2 do i=1,11

x(i)=0.08*(i-1) end do do j=1,21

y(j)=0.5+0.05*(j-1) end do do i=1,11

do j=1,21

call qiugeng(x(i),y(j),t(i,j),u(i,j)) !调用非线性方程组的牛顿迭代解法 end do end do

do i=1,6

u0(i)=0+(i-1)*0.4

t0(i)=0+(i-1)*0.2 end do

write(30,\的值为：')\

do i=1,11 do j=1,21

f(i,j)=p22(u(i,j),t(i,j),u0,t0,h,tao,z) !对于X,Y对应的U,T进行二元插值 write(*,\write(30,\

end do

end do

k=1

write(30,\收敛过程为：')\

do while(1>0)

allocate(b(11,k+1),g(21,k+1),bt(k+1,11),gt(k+1,21),btb(k+1,k+1),gtg(k+1,k+1),c(k+1,k+1)) allocate(btuu(k+1,21),uu(11,21),btbni(k+1,k+1),gtgni(k+1,k+1),bniu(k+1,21),ggni(21,k+1)) do i=1,k+1 do j=1,11

b(j,i)=x(j)**(i-1) end do end do do i=1,k+1 do j=1,21

g(j,i)=y(j)**(i-1)

end do end do

bt=.t.b !调用IMSL库函数，进行曲面拟合的矩阵计算 gt=.t.g btb=bt.x.b gtg=gt.x.g btbni=.i.btb gtgni=.i.gtg do i=1,11 do j=1,21 uu(i,j)=f(i,j) end do end do

btuu=bt.x.uu

bniu=btbni.x.btuu ggni=g.x.gtgni c=bniu.x.ggni

jingdu=0.0 do i=1,11

do j=1,21

jingdu=jingdu+(f(i,j)-pxy(x(i),y(j),c))**2 !调用多项式表达式PXY函数，计算拟合精度 end do

end do

write(*,\write(30,\if(jingdu<=0.1e-6) exit

deallocate(b,g,bt,gt,btb,gtg,c)

deallocate(btuu,uu,btbni,gtgni,bniu,ggni) k=k+1

end do

write(30,\的系数C(R,S)为：')\do i=1,k+1

write(*,\write(30,\end do

do i=1,8 xn(i)=0.1*i end do do j=1,5

yn(j)=0.5+0.2*j end do do i=1,8

do j=1,5

call qiugeng(xn(i),yn(j),t(i,j),u(i,j)) !对X*,Y*进行新的非线性方程组计算，得到新的一组T,U end do end do do i=1,8

do j=1,5

f(i,j)=p22(u(i,j),t(i,j),u0,t0,h,tao,z) !插值计算新的F(X*,Y*) end do end do

write(30,\的值为：')\do i=1,8

do j=1,5

write(*,\

e20.12,

e20.12)\

xn(i),yn(j),pxy(xn(i),yn(j),c),f(i,j)

!

输

出

X*,Y*,F(X*,Y*),P(X*,Y*)，观察逼近程度

write(30,\ end do end do

stop

end program !主程序结束

三．输出结果：

x(i),y(j),f(x,y)的值为：

0.00 0.50 0.446504069242E+00 0.00 0.55 0.324683310597E+00 0.00 0.60 0.210159730631E+00 0.00 0.65 0.103043643055E+00 0.00 0.70 0.340192524001E-02

0.00 0.75 -0.887357886680E-01 0.00 0.00 0.00 0.00

0.80 -0.173371611924E+00 0.85 -0.250534594048E+00 0.90 -0.320276491493E+00 0.95 -0.382668052743E+00

0.00 1.00 -0.437795745638E+00 0.00 1.05 -0.485758921243E+00 0.00 1.10 -0.526667236102E+00 0.00 1.15 -0.560638463024E+00 0.00 1.20 -0.587796524654E+00 0.00 1.25 -0.608269768328E+00 0.00 1.30 -0.622189446194E+00 0.00 1.35 -0.629688384281E+00 0.00 1.40 -0.630899769389E+00 0.00 1.45 -0.625956164380E+00 0.00 1.50 -0.614988550047E+00 0.08 0.50 0.638015234587E+00 0.08 0.55 0.506611795876E+00 0.08 0.60 0.382176408075E+00 0.08 0.65 0.264863525561E+00 0.08 0.70 0.154780230633E+00 0.08 0.08 0.08 0.08

0.75 0.519927097057E-01 0.80 -0.434680179481E-01 0.85 -0.131601038093E+00 0.90 -0.212431072587E+00

0.08 0.95 -0.286004537561E+00 0.08 1.00 -0.352386059629E+00 0.08 1.05 -0.411655437858E+00 0.08 1.10 -0.463904893928E+00 0.08 1.15 -0.509236708669E+00 0.08 1.20 -0.547761104059E+00 0.08 1.25 -0.579594377202E+00 0.08 1.30 -0.604857251213E+00 0.08 1.35 -0.623673426166E+00 0.08 1.40 -0.636168257103E+00 0.08 1.45 -0.642467668455E+00 0.08 1.50 -0.642697116743E+00 0.16 0.50 0.840081396883E+00 0.16 0.55 0.699764166271E+00 0.16 0.60 0.566061444836E+00 0.16 0.65 0.439171614444E+00 0.16 0.70 0.319242167288E+00 0.16 0.75 0.206376218657E+00 0.16 0.80 0.100638546919E+00

0.16 0.85 0.206075984654E-02 0.16 0.16 0.16 0.16

0.90 -0.893540080136E-01 0.95 -0.173626954709E+00 1.00 -0.250799943998E+00 1.05 -0.320932252750E+00

0.16 1.10 -0.384097719188E+00 0.16 1.15 -0.440382160815E+00 0.16 1.20 -0.489881139368E+00 0.16 1.25 -0.532697954784E+00 0.16 1.30 -0.568941871351E+00 0.16 1.35 -0.598726545062E+00 0.16 1.40 -0.622168637157E+00 0.16 1.45 -0.639386546614E+00 0.16 1.50 -0.650499364497E+00 0.24 0.50 0.105151509248E+01 0.24 0.55 0.902927427685E+00 0.24 0.60 0.760580262813E+00 0.24 0.65 0.624715195878E+00 0.24 0.70 0.495519776775E+00 0.24 0.75 0.373134067685E+00 0.24 0.80 0.257656776698E+00 0.24 0.24 0.24 0.24

0.85 0.149150588766E+00 0.90 0.476470043700E-01 0.95 -0.468493081747E-01 1.00 -0.134356747670E+00

0.24 1.05 -0.214913334079E+00 0.24 1.10 -0.288573687009E+00 0.24 1.15 -0.355406352189E+00 0.24 1.20 -0.415491385199E+00 0.24 1.25 -0.468918240380E+00 0.24 1.30 -0.515783875727E+00 0.24 1.35 -0.556191070573E+00 0.24 1.40 -0.590246929357E+00 0.24 1.45 -0.618061557834E+00 0.24 1.50 -0.639746851961E+00 0.32 0.50 0.127124675843E+01 0.32 0.55 0.111500201761E+01 0.32 0.60 0.964607722308E+00 0.32 0.65 0.820347363242E+00 0.32 0.70 0.682447673529E+00 0.32 0.75 0.551085227171E+00 0.32 0.80 0.426392408546E+00 0.32 0.85 0.308463021403E+00 0.32 0.90 0.197357157387E+00

0.32 0.95 0.931056491500E-01 0.32 0.32 0.32 0.32

1.00 -0.428596469252E-02 1.05 -0.948339139383E-01 1.10 -0.178572979137E+00 1.15 -0.255553768875E+00

0.32 1.20 -0.325840141101E+00 0.32 1.25 -0.389506980643E+00 0.32 1.30 -0.446638198547E+00 0.32 1.35 -0.497324947822E+00 0.32 1.40 -0.541664031390E+00 0.32 1.45 -0.579756487909E+00 0.32 1.50 -0.611706299458E+00 0.40 0.50 0.149832107231E+01 0.40 0.55 0.133499864071E+01 0.40 0.60 0.117712512422E+01 0.40 0.65 0.102502405047E+01 0.40 0.70 0.878960016740E+00 0.40 0.75 0.739145103470E+00 0.40 0.80 0.605744887693E+00 0.40 0.85 0.478883880984E+00 0.40 0.90 0.358650648829E+00 0.40 0.40 0.40 0.40

0.95 0.245102261159E+00 1.00 0.138268376673E+00 1.05 0.381548905914E-01 1.10 -0.552527979027E-01

0.40 1.15 -0.141986870438E+00 0.40 1.20 -0.222094430862E+00 0.40 1.25 -0.295635227198E+00 0.40 1.30 -0.362679507525E+00 0.40 1.35 -0.423306161793E+00 0.40 1.40 -0.477601035622E+00 0.40 1.45 -0.525655428548E+00 0.40 1.50 -0.567564753485E+00 0.48 0.50 0.173189277935E+01 0.48 0.55 0.156203460147E+01 0.48 0.60 0.139721693052E+01 0.48 0.65 0.123780101007E+01 0.48 0.70 0.108408753012E+01 0.48 0.75 0.936322767079E+00 0.48 0.80 0.794704459103E+00 0.48 0.85 0.659387211130E+00 0.48 0.90 0.530487603091E+00 0.48 0.95 0.408088704433E+00 0.48 1.00 0.292244222111E+00

0.48 1.05 0.182982228522E+00 0.48 0.48 0.48 0.48

1.10 0.803085151786E-01 1.15 -0.157904052160E-01 1.20 -0.105344546036E+00 1.25 -0.188398086861E+00

0.48 1.30 -0.265007147050E+00 0.48 1.35 -0.335237837830E+00 0.48 1.40 -0.399164503807E+00 0.48 1.45 -0.456868144112E+00 0.48 1.50 -0.508435001840E+00 0.56 0.50 0.197122185031E+01 0.56 0.55 0.179532964555E+01 0.56 0.60 0.162406714373E+01 0.56 0.65 0.145783060027E+01 0.56 0.70 0.129695465723E+01 0.56 0.75 0.114171810591E+01 0.56 0.80 0.992349529967E+00 0.56 0.85 0.849032670029E+00 0.56 0.90 0.711911361038E+00 0.56 0.95 0.581094170177E+00 0.56 1.00 0.456658526807E+00 0.56 0.56 0.56 0.56

1.05 1.10 1.15 1.20

0.338654511410E+00 0.227108271757E+00 0.122025104649E+00 0.233922218851E-01

0.56 1.25 -0.688187015106E-01 0.56 1.30 -0.154649345159E+00 0.56 1.35 -0.234152668459E+00 0.56 1.40 -0.307391094659E+00 0.56 1.45 -0.374434865587E+00 0.56 1.50 -0.435360560081E+00 0.64 0.50 0.221566795781E+01 0.64 0.55 0.203420120718E+01 0.64 0.60 0.185695519841E+01 0.64 0.65 0.168435820226E+01 0.64 0.70 0.151677637624E+01 0.64 0.75 0.135451905345E+01 0.64 0.80 0.119784409041E+01 0.64 0.85 0.104696305150E+01 0.64 0.90 0.902046085546E+00 0.64 0.95 0.763226480496E+00 0.64 1.00 0.630604826592E+00 0.64 1.05 0.504252821144E+00 0.64 1.10 0.384216723598E+00

0.64 1.15 0.270520485377E+00 0.64 0.64 0.64 0.64

1.20 0.163168580672E+00 1.25 0.621485545494E-01 1.30 -0.325666241909E-01 1.35 -0.121016540527E+00

0.64 1.40 -0.203251405858E+00 0.64 1.45 -0.279330366030E+00 0.64 1.50 -0.349319963960E+00 0.72 0.50 0.246468435172E+01 0.72 0.55 0.227805908576E+01 0.72 0.60 0.209525133562E+01 0.72 0.65 0.191671819269E+01 0.72 0.70 0.174285467524E+01 0.72 0.75 0.157399845779E+01 0.72 0.80 0.141043485221E+01 0.72 0.85 0.125240175695E+01 0.72 0.90 0.110009440601E+01 0.72 0.95 0.953669846066E+00 0.72 1.00 0.813251050248E+00 0.72 1.05 0.678930739188E+00 0.72 1.10 0.550774846333E+00 0.72 0.72 0.72 0.72

1.15 1.20 1.25 1.30

0.428825676240E+00 0.313104771814E+00 0.203615505489E+00 0.100345468840E+00

0.72 1.35 0.326855518060E-02 0.72 1.40 -0.876530762377E-01 0.72 1.45 -0.172467258156E+00 0.72 1.50 -0.251230230791E+00 0.80 0.50 0.271781127769E+01 0.80 0.55 0.252639964518E+01 0.80 0.60 0.233841150692E+01 0.80 0.65 0.215432947429E+01 0.80 0.70 0.197457463179E+01 0.80 0.75 0.179951063395E+01 0.80 0.80 0.162944825708E+01 0.80 0.85 0.146465006428E+01 0.80 0.90 0.130533496149E+01 0.80 0.95 0.115168260955E+01 0.80 1.00 0.100383740553E+01 0.80 1.05 0.861912322343E+00 0.80 1.10 0.725992357070E+00 0.80 1.15 0.596137699001E+00 0.80 1.20 0.472386617009E+00

0.80 1.25 0.354758086196E+00 0.80 0.80 0.80 0.80

1.30 0.243254169232E+00 1.35 0.137862207446E+00 1.40 0.385567552879E-01 1.45 -0.546986107068E-01

0.80 1.50 -0.141949673951E+00 收敛过程为：

1 0.322090912052E+01 2 0.465995586516E-02 3 0.172116950223E-03 4 0.330955599876E-05 5 0.254140048201E-07 P(X,Y)的系数C(R,S)为： 0.202122881195E+01 -0.366841611238E+01 0.709223888851E+00 0.848632669203E+00 -0.415911853259E+00 0.674349311206E-01 0.319186058746E+01 -0.740856939523E+00 -0.269763704957E+01 0.163168484787E+01 -0.484957990298E+00 0.606582355896E-01 0.257474299534E+00 0.157667343247E+01 -0.456481732139E+00 -0.884707745964E-01 0.105641397317E+00 -0.217009243867E-01 -0.271088067605E+00 -0.719886672981E+00 0.105360621876E+01 -0.783414768008E+00 0.290937635001E+00 -0.436313982104E-01 0.219696832599E+00 -0.191190784599E+00 -0.442238843798E-01 0.213523812579E+00 -0.126116723052E+00 0.234990953222E-01 -0.568656163532E-01

0.148748049781E+00

-0.148546543528E+00 0.537836373986E-01 -0.293390774839E-02 -0.133301432031E-02

x*(i),y*(j),pxy(x*,y*),f(x*,y*)的值为：

0.10 0.70 0.194730395299E+00 0.194720445694E+00 0.10 0.90 -0.183041813768E+00 -0.183037057652E+00 0.10 1.10 -0.445500015902E+00 -0.445497626192E+00 0.10 1.30 -0.597558839932E+00 -0.597566698361E+00 0.10 1.50 -0.646446110099E+00 -0.646459608525E+00 0.20 0.70 0.405989574383E+00 0.405979229893E+00 0.20 0.90 -0.225210780232E-01 -0.225159283858E-01 0.20 1.10 -0.338223995339E+00 -0.338220792377E+00 0.20 1.30 -0.544430431629E+00 -0.544437818794E+00 0.20 1.50 -0.647348013952E+00 -0.647361350834E+00 0.30 0.70 0.634787485245E+00 0.634777218201E+00 0.30 0.90 0.158796344760E+00 0.158801218632E+00 0.30 1.10 -0.207368546872E+00 -0.207365666484E+00 0.30 1.30 -0.465349901531E+00 -0.465357890239E+00 0.30 1.50 -0.620257138053E+00 -0.620270960252E+00 0.40 0.40 0.40 0.40

0.70 0.878969901892E+00 0.878960054174E+00 0.90 0.358646096786E+00 0.358650679886E+00 1.10 -0.552553941758E-01 -0.552527737602E-01 1.30 -0.362671037986E+00 -0.362679490576E+00

0.40 1.50 -0.567550573215E+00 -0.567564743834E+00 0.50 0.70 0.113662040332E+01 0.113661095866E+01 0.50 0.90 0.574975896592E+00 0.574980397067E+00 0.50 1.10 0.115989371867E+00 0.115992430481E+00 0.50 1.30 -0.238560390538E+00 -0.238568277504E+00 0.50 1.50 -0.491420882234E+00 -0.491434379464E+00 0.60 0.70 0.140605076097E+01 0.140604187520E+01 0.60 0.90 0.805937356411E+00 0.805941552321E+00 0.60 1.10 0.304425884660E+00 0.304429277291E+00 0.60 1.30 -0.950089136211E-01 -0.950160987644E-01 0.60 1.50 -0.393889814668E+00 -0.393902288570E+00 0.70 0.70 0.168579132612E+01 0.168578362972E+01 0.70 0.90 0.104987779953E+01 0.104988121574E+01 0.70 1.10 0.508291093378E+00 0.508293839537E+00 0.70 1.30 0.661563897377E-01 0.661488327126E-01 0.70 1.50 -0.276822021603E+00 -0.276834317930E+00 0.80 0.70 0.197458141737E+01 0.197457471909E+01 0.80 0.90 0.130533207763E+01 0.130533503920E+01 0.80 1.10 0.725989354460E+00 0.725992423865E+00

0.80 1.30 0.243260824222E+00 0.243254224234E+00 0.80 1.50 -0.141938770422E+00 -0.141949631258E+00

四、程序运行结果显示：

本文来源：https://www.bwwdw.com/article/laf.html