三种有限元差分法求散度

云卷云舒

各位大神，新手求助求解散度！有经向纬向两组速度数据，（附加图已给出,这是我读出数据后又整理了一下输出的，在我读取源文件是是先读的一行x）
用前差、后差、中心差求解各点的散度（偏u偏x+偏v偏y）（由于边界没法求，有的点没有散度）
但是我求出来三种方法的散度值完全不一样，不知道问题出在那里，求高手指点一下，万分感激！！！
一下是我的程序：（只看计算部分即可，print部分可以不看）

program test11

   implicit none

   integer::n,i,j

   real::x(17),y(17)

   real,parameter::d=0.25

   real::u(17,17),v(17,17)

   real::df_ux(16,16),df_vy(16,16),df(16,16)

   real::da_ux(2:17,2:17),da_vy(2:17,2:17),da(2:17,2:17)

   real::dm_ux(2:16,2:16),dm_vy(2:16,2:16),dm(2:16,2:16)!----------------一种方法不能只开一个数组，因为两次求导的内外层循环顺序反了！！！！

   open(unit=10,file='u.txt',status='old')

   open(unit=20,file='v.txt',status='old')

   read(10,*) (x(i),i=1,17)  !-----------注意这里

   do i=1,17

       read(20,'(18f12.7)') y(i),v(i,:)

       read(10,'(18f12.7)') y(i),u(i,:)      

   end do

   !============================计算====================================

   !=============================前差=======================

   !---------u/x

    do j=1,16

       do i=1,16

           df_ux(i,j)=(u(i+1,j)-u(i,j))/d!----i对应x,j对应y

       end do

   end do

   !---------v/y

   do i=1,16

       do j=1,16

           df_vy(i,j)=(v(i,j+1)-v(i,j))/d

       end do

   end do

     df=df_ux+df_vy

   !=============================后差=======================

    !---------u/x

   do j=2,17

       do i=17,2,-1

           da_ux(i,j)=(u(i-1,j)-u(i,j))/d

       end do

   end do

   !---------v/y

   do i=2,17

       do j=17,2,-1

           da_vy(i,j)=(v(i,j-1)-v(i,j))/d

       end do

   end do

     da=da_ux+da_vy

   !===========================中心差=======================

   !---------u/x

   do j=2,16

       do i=2,16

           dm_ux(i,j)=(u(i+1,j)-u(i-1,j))/(2.0*d)

       end do

   end do

   !---------v/y

   do j=2,16

       do i=2,16

           dm_vy(i,j)=(v(i,j+1)-v(i,j-1))/(2.0*d)

       end do

   end do

     dm=dm_ux+dm_vy

   !==========================print===========================================

   !=============================前差=======================

   open(unit=11,file='div_f.txt')

   do i=16,1,-1

       write(11,'(18f12.7)') y(i),df(i,:)  

   end do

   write(11,'(12x,16f12.7)') (x(i),i=1,16)

   !-------------------------------------

   print*,"前差："

   do i=16,1,-1

       print'(18f10.7)', y(i),df(i,:)   

   end do

   print'(10x,16f10.7)',(x(i),i=1,16)

   !=============================后差=======================

   open(unit=12,file='div_a.txt')

   do i=17,2,-1

       write(12,'(18f12.7)')y(i),da(i,:)   

   end do

   write(12,'(12x,16f12.7)')(x(i),i=2,17)

   !---------------------------------------

     print*,"后差："

   do i=17,2,-1

       print'(18f10.7)', y(i),da(i,:)   

   end do

   print'(10x,16f10.7)',(x(i),i=2,17)

   !===========================中心差=======================

   open(unit=13,file='div_m.txt')

   do i=16,2,-1

       write(13,'(18f12.7)') y(i),dm(i,:)  

   end do

   write(13,'(12x,15f12.7)') (x(i),i=2,16)

    !-------------------------------

     print*,"中心差："

   do i=16,2,-1

       print'(18f10.7)', y(i),dm(i,:)   

   end do

   print'(10x,15f10.7)', (x(i),i=2,16)

   end program

li913

各种差分格式的结果有差异，这是肯定的，不然就没必要那么多格式了。
后差写错了，符号不对。