您所在位置：主题：程序无误，输出结果不对??????????
wangxinling
[专家分：0] 发布于 2012-07-23 19:26:00
模型：三层介质，速度分别是1600m/s,3300m/s,5000m/s，800*800网格，网格间距1m*1m。（网格划分就是知道每个节点v(I,j)上的速度。）程序目的是：求任意一点（不止节点）的速度和导数vv,vvx,vvz。公式g(x,z)=1/4*(dz)^cdx，其中z,d,c,d,x都是矩阵。
dimension v(1000,1000),vx(1000,1000),vz(1000,1000)
*,d(4,4),dt(4,4),dd(4,4),ddt(4,4),c(1000,1000)
*,x(1000,1000),z(1000,1000),zt(1000,1000)
*,xx(1000,1000),zz(1000,1000),zzt(1000,1000)
*,ztdt(1000,1000),zztddt(1000,1000),dx(1000,1000),ddxx(1000,1000)
*,ztdtc(1000,1000),zztddtc(1000,1000)
*,ztdtcdx(1,1),ztdtcddxx(1,1),zztddtcdx(1,1)
real v,vx,vz,d,dt,c,x,z,zt,xx,zz,zzt,vv,vvx,vvz,dx,ddx,ddz,h,l
*,xxx,zzz,xxxx,zzzz
h=1
L=1
do i=1,801
do j=1,201
v(i,j)=1600
enddo
enddo
do i=1,801
do j=202,401
v(i,j)=3300
enddo
enddo
do i=1,801
do j=402,801
v(i,j)=5000
enddo
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc输出所有节点上的速度
do i=1,801
do j=1,801
c print*,i,j,v(i,j)
enddo
enddo
c 内部节点导数
do i=2,800
do j=2,800
vx(i,j)=(v(i+1,j)-v(i-1,j))/2*h
vz(i,j)=(v(i,j+1)-v(i,j-1))/2*l
enddo
enddo
c 边界节点导数
do j=1,800
vx(1,j)=(v(2,j)-v(1,j))/h
vz(1,j)=(v(1,j+1)-v(1,j))/l
vx(801,j)=(v(801,j)-v(800,j))/h
vz(801,j)=(v(801,j+1)-v(801,j))/l
enddo
do i=1,800
vx(i,1)=(v(i+1,1)-v(i,1))/h
vz(i,1)=(v(i,2)-v(i,1))/l
vx(i,801)=(v(i+1,801)-v(i,801))/h
vz(i,801)=(v(i,801)-v(i,800))/l
enddo
vx(801,801)=(v(801,801)-v(800,801))/h
vz(801,801)=(v(801,801)-v(801,800))/l
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc输出所有节点上的导数
do i=1,801
do j=1,801
c print*,i,j,vx(i,j),vz(i,j)
enddo
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc插值：求任意一点的速度及导数
read*,xxxx,zzzz
if(xxxx.ge.1.and.xxxx.le.798.and.zzzz.ge.1.and.zzzz.le.798)then
i=int(xxxx/h)+1
j=int(zzzz/l)+1
xxx=(i-1)*h
zzz=(j-1)*l
if(xxxx.gt.xxx.and.xxxx.lt.xxx+h.and.zzzz.gt.zzz.and.zzzz.lt.zzz+L)then
d(1,1)=-1
d(1,2)=2
d(1,3)=-1
d(1,4)=0
d(2,1)=3
d(2,2)=-5
d(2,3)=0
d(2,4)=2
d(3,1)=-3
d(3,2)=4
d(3,3)=1
d(3,4)=0
d(4,1)=1
d(4,2)=-1
d(4,3)=0
d(4,4)=0
dt(1,1)=-1
dt(1,2)=3
dt(1,3)=-3
dt(1,4)=1
dt(2,1)=2
dt(2,2)=-5
dt(2,3)=4
dt(2,4)=-1
dt(3,1)=-1
dt(3,2)=0
dt(3,3)=1
dt(3,4)=0
dt(4,1)=0
dt(4,2)=2
dt(4,3)=0
dt(4,4)=0
dd(1,1)=0
dd(1,2)=-3
dd(1,3)=4
dd(1,4)=-1
dd(2,1)=0
dd(2,2)=9
dd(2,3)=-10
dd(2,4)=0
dd(3,1)=0
dd(3,2)=-9
dd(3,3)=8
dd(3,4)=1
dd(4,1)=0
dd(4,2)=3
dd(4,3)=-2
dd(4,4)=0
ddt(1,1)=0
ddt(1,2)=0
ddt(1,3)=0
ddt(1,4)=0
ddt(2,1)=-3
ddt(2,2)=9
ddt(2,3)=-9
ddt(2,4)=3
ddt(3,1)=4
ddt(3,2)=-10
ddt(3,3)=8
ddt(3,4)=-2
ddt(4,1)=-1
ddt(4,2)=0
ddt(4,3)=1
ddt(4,4)=0
c(1,1)=v(i-1,j-1)
c(1,2)=v(i,j-1)
c(1,3)=v(i+1,j-1)
c(1,4)=v(i+2,j-1)
c(2,1)=v(i-1,j)
c(2,2)=v(i,j)
c(2,3)=v(i+1,j)
c(2,4)=v(i+2,j)
c(3,1)=v(i-1,j+1)
c(3,2)=v(i,j+1)
c(3,3)=v(i+1,j+1)
c(3,4)=v(i+2,j+1)
c(4,1)=v(i-1,j+2)
c(4,2)=v(i,j+2)
c(4,3)=v(i+1,j+2)
c(4,4)=v(i+2,j+2)
ddx=(xxxx-xxx)/h
ddz=(zzzz-zzz)/l
x(1,1)=ddx**3
x(2,1)=ddx**2
x(3,1)=ddx**1
x(4,1)=ddx**0
z(1,1)=ddz**3
z(2,1)=ddz**2
z(3,1)=ddz**1
z(4,1)=ddz**0
zt(1,1)=z(1,1)
zt(1,2)=z(2,1)
zt(1,3)=z(3,1)
zt(1,4)=z(4,1)
xx(1,1)=ddx**3
xx(2,1)=ddx**2
xx(3,1)=ddx**1
xx(4,1)=ddx**0
zz(1,1)=ddz**3
zz(2,1)=ddz**2
zz(3,1)=ddz**1
zz(4,1)=ddz**0
zzt(1,1)=zz(1,1)
zzt(1,2)=zz(2,1)
zzt(1,3)=zz(3,1)
zzt(1,4)=zz(4,1)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 两个矩阵相乘公式，ii为行，jj为列。
c ztdt=zt*dt
do ii=1,1
do jj=1,4
ztdt(ii,jj)=0
do k=1,4
ztdt(ii,jj)=ztdt(ii,jj)+zt(ii,k)*dt(k,jj)
enddo
enddo
enddo
c zztddt=zzt*ddt
do ii=1,1
do jj=1,4
zztddt(ii,jj)=0
do k=1,4
zztddt(ii,jj)=zztddt(ii,jj)+zzt(ii,k)*ddt(k,jj)
enddo
enddo
enddo
c dx=d*x
do ii=1,4
do jj=1,1
dx(ii,jj)=0
do k=1,4
dx(ii,jj)=dx(ii,jj)+d(ii,k)*x(k,jj)
enddo
enddo
enddo
c ddxx=dd*xx
do ii=1,4
do jj=1,1
ddxx(ii,jj)=0
do k=1,4
ddxx(ii,jj)=ddxx(ii,jj)+dd(ii,k)*xx(k,jj)
enddo
enddo
enddo
c ztdtc=zt*dt*c=ztdt*c
do ii=1,1
do jj=1,4
ztdtc(ii,jj)=0
do k=1,4
ztdtc(ii,jj)=ztdtc(ii,jj)+ztdt(ii,k)*c(k,jj)
enddo
enddo
enddo
c zztddtc=zzt*ddt*c=zztddt*c
do ii=1,1
do jj=1,4
zztddtc(ii,jj)=0
do k=1,4
zztddtc(ii,jj)=zztddtc(ii,jj)+zztddt(ii,k)*c(k,jj)
enddo
enddo
enddo
c ztdtcdx=zt*dt*c*dx=ztdtc*dx
do ii=1,1
do jj=1,1
ztdtcdx(ii,jj)=0
do k=1,4
ztdtcdx(ii,jj)=ztdtcdx(ii,jj)+ztdtc(ii,k)*dx(k,jj)
vv=1/4*ztdtcdx(ii,jj)
enddo
enddo
enddo
c ztdtcddxx=zt*dt*c*dd*xx=ztdtc*ddxx
do ii=1,1
do jj=1,1
ztdtcddxx(ii,jj)=0
do k=1,4
ztdtcddxx(ii,jj)=ztdtcddxx(ii,jj)+ztdtc(ii,k)*ddxx(k,jj)
vvx=1/4*ztdtcddxx(ii,jj)
enddo
enddo
enddo
c zztddtcdx=zzt*ddt*c*d*x=zztddtc*dx
do ii=1,1
do jj=1,1
zztddtcdx(ii,jj)=0
do k=1,4
zztddtcdx(ii,jj)=zztddtcdx(ii,jj)+zztddtc(ii,k)*dx(k,jj)
vvz=1/4*zztddtcdx(ii,jj)
enddo
enddo
enddo
elseif(zzzz.eq.zzz.and.xxxx.ge.xxx.and.xxxx.le.xxx+h)then
vv=(xxxx-xxx)/h*v(i+1,j)+(xxx+h-xxxx)/h*v(i,j)
vvx=(xxxx-xxx)/h*vx(i+1,j)+(xxx+h-xxxx)/h*vx(i,j)
vvz=(xxxx-xxx)/h*vz(i+1,j)+(xxx+h-xxxx)/h*vz(i,j)
elseif(zzzz.eq.zzz+l.and.xxxx.ge.xxx.and.xxxx.le.xxx+h)then
vv=(xxxx-xxx)/h*v(i+1,j+1)+(xxx+h-xxxx)/h*v(i,j+1)
vvx=(xxxx-xxx)/h*vx(i+1,j+1)+(xxx+h-xxxx)/h*vx(i,j+1)
vvz=(xxxx-xxx)/h*vz(i+1,j+1)+(xxx+h-xxxx)/h*vz(i,j+1)
elseif(xxxx.eq.xxx.and.zzzz.ge.zzz.and.zzzz.le.zzz+l)then
vv=(zzzz-zzz)/l*v(i,j+1)+(zzz+L-zzzz)/l*v(i,j)
vvx=(zzzz-zzz)/l*vx(i,j+1)+(zzz+L-zzzz)/l*vx(i,j)
vvz=(zzzz-zzz)/l*vz(i,j+1)+(zzz+L-zzzz)/l*vz(i,j)
elseif(xxxx.eq.xxx+h.and.zzzz.ge.zzz.and.zzzz.le.zzz+l)then
vv=(zzzz-zzz)/l*v(i+1,j+1)+(zzz+L-zzzz)/l*v(i+1,j)
vvx=(zzzz-zzz)/l*vx(i+1,j+1)+(zzz+L-zzzz)/l*vx(i+1,j)
vvz=(zzzz-zzz)/l*vz(i+1,j+1)+(zzz+L-zzzz)/l*vz(i+1,j)
endif
else
i=int(xxxx/h)+1
j=int(zzzz/l)+1
xxx=(i-1)*h
zzz=(j-1)*l
v1=(xxxx-xxx)/h*v(i+1,j)+(xxx+h-xxxx)/h*v(i,j)
v2=(xxxx-xxx)/h*v(i+1,j+1)+(xxx+h-xxxx)/h*v(i,j+1)
vv=(zzzz-zzz)/L*v2+(zzz+l-zzzz)/L*v1
v1x=(xxxx-xxx)/h*vx(i+1,j)+(xxx+h-xxxx)/h*vx(i,j)
v1z=(xxxx-xxx)/h*vz(i+1,j)+(xxx+h-xxxx)/h*vz(i,j)
v2x=(xxxx-xxx)/h*vx(i+1,j+1)+(xxx+h-xxxx)/h*vx(i,j+1)
v2z=(xxxx-xxx)/h*vz(i+1,j+1)+(xxx+h-xxxx)/h*vz(i,j+1)
vvx=(zzzz-zzz)/L*v2x+(zzz+l-zzzz)/L*v1x
vvz=(zzzz-zzz)/L*v2z+(zzz+l-zzzz)/L*v1z
endif
print*,vv,vvx,vvz
end
dimension v(1000,1000),vx(1000,1000),vz(1000,1000)
*,d(4,4),dt(4,4),dd(4,4),ddt(4,4),c(1000,1000)
*,x(1000,1000),z(1000,1000),zt(1000,1000)
*,xx(1000,1000),zz(1000,1000),zzt(1000,1000)
*,ztdt(1000,1000),zztddt(1000,1000),dx(1000,1000),ddxx(1000,1000)
*,ztdtc(1000,1000),zztddtc(1000,1000)
*,ztdtcdx(1,1),ztdtcddxx(1,1),zztddtcdx(1,1)
real v,vx,vz,d,dt,c,x,z,zt,xx,zz,zzt,vv,vvx,vvz,dx,ddx,ddz,h,l
*,xxx,zzz,xxxx,zzzz
h=1
L=1
do i=1,801
do j=1,201
v(i,j)=1600
enddo
enddo
do i=1,801
do j=202,401
v(i,j)=3300
enddo
enddo
do i=1,801
do j=402,801
v(i,j)=5000
enddo
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc输出所有节点上的速度
do i=1,801
do j=1,801
c print*,i,j,v(i,j)
enddo
enddo
c 内部节点导数
do i=2,800
do j=2,800
vx(i,j)=(v(i+1,j)-v(i-1,j))/2*h
vz(i,j)=(v(i,j+1)-v(i,j-1))/2*l
enddo
enddo
c 边界节点导数
do j=1,800
vx(1,j)=(v(2,j)-v(1,j))/h
vz(1,j)=(v(1,j+1)-v(1,j))/l
vx(801,j)=(v(801,j)-v(800,j))/h
vz(801,j)=(v(801,j+1)-v(801,j))/l
enddo
do i=1,800
vx(i,1)=(v(i+1,1)-v(i,1))/h
vz(i,1)=(v(i,2)-v(i,1))/l
vx(i,801)=(v(i+1,801)-v(i,801))/h
vz(i,801)=(v(i,801)-v(i,800))/l
enddo
vx(801,801)=(v(801,801)-v(800,801))/h
vz(801,801)=(v(801,801)-v(801,800))/l
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc输出所有节点上的导数
do i=1,801
do j=1,801
c print*,i,j,vx(i,j),vz(i,j)
enddo
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc插值：求任意一点的速度及导数
read*,xxxx,zzzz
if(xxxx.ge.1.and.xxxx.le.798.and.zzzz.ge.1.and.zzzz.le.798)then
i=int(xxxx/h)+1
j=int(zzzz/l)+1
xxx=(i-1)*h
zzz=(j-1)*l
if(xxxx.gt.xxx.and.xxxx.lt.xxx+h.and.zzzz.gt.zzz.and.zzzz.lt.zzz+L)then
d(1,1)=-1
d(1,2)=2
d(1,3)=-1
d(1,4)=0
d(2,1)=3
d(2,2)=-5
d(2,3)=0
d(2,4)=2
d(3,1)=-3
d(3,2)=4
d(3,3)=1
d(3,4)=0
d(4,1)=1
d(4,2)=-1
d(4,3)=0
d(4,4)=0
dt(1,1)=-1
dt(1,2)=3
dt(1,3)=-3
dt(1,4)=1
dt(2,1)=2
dt(2,2)=-5
dt(2,3)=4
dt(2,4)=-1
dt(3,1)=-1
dt(3,2)=0
dt(3,3)=1
dt(3,4)=0
dt(4,1)=0
dt(4,2)=2
dt(4,3)=0
dt(4,4)=0
dd(1,1)=0
dd(1,2)=-3
dd(1,3)=4
dd(1,4)=-1
dd(2,1)=0
dd(2,2)=9
dd(2,3)=-10
dd(2,4)=0
dd(3,1)=0
dd(3,2)=-9
dd(3,3)=8
dd(3,4)=1
dd(4,1)=0
dd(4,2)=3
dd(4,3)=-2
dd(4,4)=0
ddt(1,1)=0
ddt(1,2)=0
ddt(1,3)=0
ddt(1,4)=0
ddt(2,1)=-3
ddt(2,2)=9
ddt(2,3)=-9
ddt(2,4)=3
ddt(3,1)=4
ddt(3,2)=-10
ddt(3,3)=8
ddt(3,4)=-2
ddt(4,1)=-1
ddt(4,2)=0
ddt(4,3)=1
ddt(4,4)=0
c(1,1)=v(i-1,j-1)
c(1,2)=v(i,j-1)
c(1,3)=v(i+1,j-1)
c(1,4)=v(i+2,j-1)
c(2,1)=v(i-1,j)
c(2,2)=v(i,j)
c(2,3)=v(i+1,j)
c(2,4)=v(i+2,j)
c(3,1)=v(i-1,j+1)
c(3,2)=v(i,j+1)
c(3,3)=v(i+1,j+1)
c(3,4)=v(i+2,j+1)
c(4,1)=v(i-1,j+2)
c(4,2)=v(i,j+2)
c(4,3)=v(i+1,j+2)
c(4,4)=v(i+2,j+2)
ddx=(xxxx-xxx)/h
ddz=(zzzz-zzz)/l
x(1,1)=ddx**3
x(2,1)=ddx**2
x(3,1)=ddx**1
x(4,1)=ddx**0
z(1,1)=ddz**3
z(2,1)=ddz**2
z(3,1)=ddz**1
z(4,1)=ddz**0
zt(1,1)=z(1,1)
zt(1,2)=z(2,1)
zt(1,3)=z(3,1)
zt(1,4)=z(4,1)
xx(1,1)=ddx**3
xx(2,1)=ddx**2
xx(3,1)=ddx**1
xx(4,1)=ddx**0
zz(1,1)=ddz**3
zz(2,1)=ddz**2
zz(3,1)=ddz**1
zz(4,1)=ddz**0
zzt(1,1)=zz(1,1)
zzt(1,2)=zz(2,1)
zzt(1,3)=zz(3,1)
zzt(1,4)=zz(4,1)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 两个矩阵相乘公式，ii为行，jj为列。
c ztdt=zt*dt
do ii=1,1
do jj=1,4
ztdt(ii,jj)=0
do k=1,4
ztdt(ii,jj)=ztdt(ii,jj)+zt(ii,k)*dt(k,jj)
enddo
enddo
enddo
c zztddt=zzt*ddt
do ii=1,1
do jj=1,4
zztddt(ii,jj)=0
do k=1,4
zztddt(ii,jj)=zztddt(ii,jj)+zzt(ii,k)*ddt(k,jj)
enddo
enddo
enddo
c dx=d*x
do ii=1,4
do jj=1,1
dx(ii,jj)=0
do k=1,4
dx(ii,jj)=dx(ii,jj)+d(ii,k)*x(k,jj)
enddo
enddo
enddo
c ddxx=dd*xx
do ii=1,4
do jj=1,1
ddxx(ii,jj)=0
do k=1,4
ddxx(ii,jj)=ddxx(ii,jj)+dd(ii,k)*xx(k,jj)
enddo
enddo
enddo
c ztdtc=zt*dt*c=ztdt*c
do ii=1,1
do jj=1,4
ztdtc(ii,jj)=0
do k=1,4
ztdtc(ii,jj)=ztdtc(ii,jj)+ztdt(ii,k)*c(k,jj)
enddo
enddo
enddo
c zztddtc=zzt*ddt*c=zztddt*c
do ii=1,1
do jj=1,4
zztddtc(ii,jj)=0
do k=1,4
zztddtc(ii,jj)=zztddtc(ii,jj)+zztddt(ii,k)*c(k,jj)
enddo
enddo
enddo
c ztdtcdx=zt*dt*c*dx=ztdtc*dx
do ii=1,1
do jj=1,1
ztdtcdx(ii,jj)=0
do k=1,4
ztdtcdx(ii,jj)=ztdtcdx(ii,jj)+ztdtc(ii,k)*dx(k,jj)
vv=1/4*ztdtcdx(ii,jj)
enddo
enddo
enddo
c ztdtcddxx=zt*dt*c*dd*xx=ztdtc*ddxx
do ii=1,1
do jj=1,1
ztdtcddxx(ii,jj)=0
do k=1,4
ztdtcddxx(ii,jj)=ztdtcddxx(ii,jj)+ztdtc(ii,k)*ddxx(k,jj)
vvx=1/4*ztdtcddxx(ii,jj)
enddo
enddo
enddo
c zztddtcdx=zzt*ddt*c*d*x=zztddtc*dx
do ii=1,1
do jj=1,1
zztddtcdx(ii,jj)=0
do k=1,4
zztddtcdx(ii,jj)=zztddtcdx(ii,jj)+zztddtc(ii,k)*dx(k,jj)
vvz=1/4*zztddtcdx(ii,jj)
enddo
enddo
enddo
elseif(zzzz.eq.zzz.and.xxxx.ge.xxx.and.xxxx.le.xxx+h)then
vv=(xxxx-xxx)/h*v(i+1,j)+(xxx+h-xxxx)/h*v(i,j)
vvx=(xxxx-xxx)/h*vx(i+1,j)+(xxx+h-xxxx)/h*vx(i,j)
vvz=(xxxx-xxx)/h*vz(i+1,j)+(xxx+h-xxxx)/h*vz(i,j)
elseif(zzzz.eq.zzz+l.and.xxxx.ge.xxx.and.xxxx.le.xxx+h)then
vv=(xxxx-xxx)/h*v(i+1,j+1)+(xxx+h-xxxx)/h*v(i,j+1)
vvx=(xxxx-xxx)/h*vx(i+1,j+1)+(xxx+h-xxxx)/h*vx(i,j+1)
vvz=(xxxx-xxx)/h*vz(i+1,j+1)+(xxx+h-xxxx)/h*vz(i,j+1)
elseif(xxxx.eq.xxx.and.zzzz.ge.zzz.and.zzzz.le.zzz+l)then
vv=(zzzz-zzz)/l*v(i,j+1)+(zzz+L-zzzz)/l*v(i,j)
vvx=(zzzz-zzz)/l*vx(i,j+1)+(zzz+L-zzzz)/l*vx(i,j)
vvz=(zzzz-zzz)/l*vz(i,j+1)+(zzz+L-zzzz)/l*vz(i,j)
elseif(xxxx.eq.xxx+h.and.zzzz.ge.zzz.and.zzzz.le.zzz+l)then
vv=(zzzz-zzz)/l*v(i+1,j+1)+(zzz+L-zzzz)/l*v(i+1,j)
vvx=(zzzz-zzz)/l*vx(i+1,j+1)+(zzz+L-zzzz)/l*vx(i+1,j)
vvz=(zzzz-zzz)/l*vz(i+1,j+1)+(zzz+L-zzzz)/l*vz(i+1,j)
endif
else
i=int(xxxx/h)+1
j=int(zzzz/l)+1
xxx=(i-1)*h
zzz=(j-1)*l
v1=(xxxx-xxx)/h*v(i+1,j)+(xxx+h-xxxx)/h*v(i,j)
v2=(xxxx-xxx)/h*v(i+1,j+1)+(xxx+h-xxxx)/h*v(i,j+1)
vv=(zzzz-zzz)/L*v2+(zzz+l-zzzz)/L*v1
v1x=(xxxx-xxx)/h*vx(i+1,j)+(xxx+h-xxxx)/h*vx(i,j)
v1z=(xxxx-xxx)/h*vz(i+1,j)+(xxx+h-xxxx)/h*vz(i,j)
v2x=(xxxx-xxx)/h*vx(i+1,j+1)+(xxx+h-xxxx)/h*vx(i,j+1)
v2z=(xxxx-xxx)/h*vz(i+1,j+1)+(xxx+h-xxxx)/h*vz(i,j+1)
vvx=(zzzz-zzz)/L*v2x+(zzz+l-zzzz)/L*v1x
vvz=(zzzz-zzz)/L*v2z+(zzz+l-zzzz)/L*v1z
endif
print*,vv,vvx,vvz
end