fortran 一维和二维插值.docx

《fortran 一维和二维插值.docx》由会员分享，可在线阅读，更多相关《fortran 一维和二维插值.docx（6页珍藏版）》请在三一办公上搜索。

1、fortran 一维和二维插值 program Main c 声明变量 integer t character*32 nodes real Inv,Inv1,u,v real,allocatable:X(:,:,:),Y(:,:,:),Z(:,:,:),W(:,:,:), $X0(:,:,:),Y0(:,:,:),Z0(:,:,:),W0(:,:,:),x1(:),y1(:), $w1(:),z1(:,:) c 打开已有文件11读取数据行数 nodes=F1.DAT open(11,file=nodes,status=old,form=formatted) n=0 do while(.not.

2、eof(11) read(11,*,end=102) n=n+1 102 end do write(*,*) n c 返回首地址 rewind(11) c 判断三维数组大小 do i=1,(n/2) if(i*i*i.eq.n)then t=i end if end do c 调用定义的动态三维数组并读入值 allocate(x(1:t,1:t,1:t),y(1:t,1:t,1:t),z(1:t,1:t,1:t), $w(1:t,1:t,1:t),w0(1:(t-1),1:(t-1),1:(t-1), $x0(1:(t-1),1:(t-1),1:(t-1),y0(1:(t-1),1:(t-1)

3、,1:(t-1), $z0(1:(t-1),1:(t-1),1:(t-1) do 30 i=1,t,1 do 40 j=1,t,1 do 50 k=1,t,1 read(11,100) x(k,j,i),y(k,j,i),z(k,j,i),w(k,j,i) 50 continue 40 continue 30 continue rewind(11) close(11) c 计算flac 3d平均节点和坐标值 do 65 i=1,t-1,1 do 75 j=1,t-1,1 do 85 k=1,t-1,1 w0(k,j,i)=(w(i,j,k)+w(i,j,k+1)+w(i,j+1,k)+w(i,

4、j+1,k+1)+ $w(i+1,j,k)+w(i+1,j,k+1)+w(i+1,j+1,k)+w(i+1,j+1,k+1)/8 x0(k,j,i)=(x(i,j,k)+x(i,j,k+1)+x(i,j+1,k)+x(i,j+1,k+1)+ $x(i+1,j,k)+x(i+1,j,k+1)+x(i+1,j+1,k)+x(i+1,j+1,k+1)/8 y0(k,j,i)=(y(i,j,k)+y(i,j,k+1)+y(i,j+1,k)+y(i,j+1,k+1)+ $y(i+1,j,k)+y(i+1,j,k+1)+y(i+1,j+1,k)+y(i+1,j+1,k+1)/8 z0(k,j,i)=(z(

5、i,j,k)+z(i,j,k+1)+z(i,j+1,k)+z(i,j+1,k+1)+ $z(i+1,j,k)+z(i+1,j,k+1)+z(i+1,j+1,k)+z(i+1,j+1,k+1)/8 85 continue 75 continue 65 continue c 把求解的节点值输入到文件2 OPEN (UNIT=2,FILE=F2.DAT,STATUS=NEW,ACCESS=SEQUENTIAL, $FORM=FORMATTED) do 60 i=1,t-1,1 do 70 j=1,t-1,1 do 80 k=1,t-1,1 write(2,100) x0(k,j,i),y0(k,j,

6、i),z0(k,j,i),w0(k,j,i) 80 continue 70 continue 60 continue close(2) c 调用计算一维全区间插值并输出数据 open(11,file=nodes,status=old,form=formatted) allocate (x1(1:n),y1(1:n),w1(1:n) do i=1,n,1 read(11,100) x1(i),y1(i),w1(i) end do close(11) c 把数据输入到文件3 OPEN(UNIT=3,FILE=F3.DAT,STATUS=NEW,ACCESS=SEQUENTIAL, $FORM=FO

7、RMATTED) do i=1,n,1 call flac11(x1,y1,n,w1(i),InV) write(3,(e26.6) InV end do CLOSE(3) c 调用计算一维三点插值并输出数据 c 把数据输入到文件4 OPEN(UNIT=4,FILE=F4.DAT,STATUS=NEW,ACCESS=SEQUENTIAL, $FORM=FORMATTED) do i=1,n,1 call flac13(x1,y1,n,w1(i),InV1) write(4,300) x1(i),y1(i),Inv1 end do CLOSE(4) 100 format(e12.6,3e16.6

8、) 300 format(1x,x=,e12.6,3x,y=,e12.6,3x,z=,e12.6) c 调用计算二维三点插值并输出数据 c 把数据输入到文件5 allocate (z1(1:n,1:n) do j=1,n,1 do i=1,n,1 z1(i,j)=exp(-(x1(i)-y1(j) end do end do u=0.9 v=0.8 call flac2(x1,y1,z1,n,n,u,v,w) write(*,41) u,v,w 41 format(1x,x=,f7.3,1x,y=,f7.3,1x,z(x,y)=,D15.6) u=0.3 v=0.9 call flac2(x1

9、,y1,z1,n,n,u,v,w) write(*,41) u,v,w end c 调用一维全区间插值的子程序 subroutine flac11(X,Y,n,InsertX,InsertV) dimension x(n),y(n) real InsertX,InsertV,temp,X,Y InsertV=0.0 do 20 i=1,n temp=1.0 do 10 j=1,n if (i.NE.j) then temp=temp*(InsertX-X(j)/(x(i)-X(j) end if 10 continue InsertV=InsertV+y(i)*temp 20 continue

10、 end subroutine c 调用一维三点插值的子程序 subroutine flac13(x,y,n,t,z) dimension x(n),y(n) real x,y,t,z,s z=0.0 if(n.le.0)return if(n.eq.1)then z=y(1) return end if if(n.eq.2)then z=(y(1)*(t-x(2)-y(2)*(t-x(1)/(x(1)-x(2) return end if if(t.le.x(2)then k=1 m=3 else if(t.ge.x(n-1)then k=n-2 m=n else k=1 m=n 13 if

11、(iabs(k-m).ne.1)then L=(k+m)/2 if(t.lt.x(L)then m=L else k=L end if goto 13 end if if(abs(t-x(k).lt.abs(t-x(m)then k=k-1 else m=m+1 end if end if z=0.0 do 33 i=k,m s=1.0 do 23 j=k,m if (j.ne.i)then s=s*(t-x(j)/(x(i)-x(j) end if 23 continue z=z+s*y(i) 33 continue return end subroutine c 调用二维三点插值子程序 s

12、ubroutine flac2(x,y,z,n,m,u,v,w) dimension x(n),y(m),z(n,m),b(3) reaL x,y,z,u,v,w,b,hh nn=3 if(n.le.3)then ip=1 nn=n eLse if(u.le.x(2)then ip=1 eLse if(u.ge.x(n-1)then ip=n-2 eLse i=1 j=n 12 if(iabs(i-j).ne.1)then L=(i+j)/2 if(u.lt.x(L)then j=L eLse i=L end if goto 12 end if if(abs(u-x(i).lt.abs(u-x

13、(j)then ip=i-1 eLse ip=i end if end if mm=3 if(m.Le.3)then iq=1 mm=m eLse if(v.Le.y(2)then iq=1 eLse if(v.ge.y(m-1)then iq=m-2 eLse i=1 j=m 22 if(iabs(i-j).ne.1)then L=(i+j)/2 if(v.lt.y(L)then j=L eLse i=L end if goto 22 end if if(abs(v-y(i).lt.abs(v-y(j)then iq=i-1 eLse iq=i end if end if do 52 i=1

14、,nn b(i)=0.0 do 42 j=1,mm hh=z(ip+i-1,iq+j-1) do 32 k=1,mm if(k.ne.j)then hh=hh*(v-y(iq+k-1)/(y(iq+j-1)-y(iq+k-1) end if 32 continue b(i)=b(i)+hh 42 continue 52 continue w=0.0 do 72 i=1,nn hh=b(i) do 62 j=1,nn if(i.ne.j)then hh=hh*(u-x(ip+j-1)/(x(ip+i-1)-x(ip+j-1) end if 62 continue w=w+hh 72 continue return end subroutine