北航数值分析大作业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

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/l2dv.html