《数值分析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