王培杰 发表于 2014-4-20 19:53:53

MT一维正演，用SOR法不收敛

原理如图，结果一直不收敛
代码如下program main
implicit none
integer,parameter::n=1000                     !地上n层 地下n层
complex*16,parameter::i=(0d0,1d0)               
real*16,parameter::PI=3.1415926535897932384626
real*16,parameter::miu=4d0*PI*1E-7            
real*16,parameter::w=1.5                        !SOR迭代因子
integer::j,k
real*16::o(-n-1:n),o1(-n-1:n),dz(-n-1:n+1),s,f,E14,ro,fa
!o地层电导率，o1分界面电导率，dz层厚度，f频率，E14四分之一厚度处电场，ro视电阻率，fa相位
complex*16::E(-n-1:n+1),H(-n-1:n),E1(-n-1:n+1),Hs,a,b
!E电场，H磁场，Hs地表磁场
open(123,file='a.dat')
dz=1d0                  !单层厚度1m
dz(-n-1)=1.0E6            !最上层厚度1000km，视作无限远
dz(n)=1.0E6               !最底层厚度1000km，视作无限远
forall (k=-n-1:-1:1) o(k)=1E-10      !空气中电导率1E-10，视作无限大
forall(k=0:300:1) o(k)=0.1d0         !0到300米，电导率0.1
forall(k=301:600:1) o(k)=0.005d0   !301到600米，电导率0.005
forall(k=601:n:1) o(k)=0.1d0         !601到无穷远，电导率0.1
do k=-n,n
    o1(k)=(o(k)*dz(k)+o(k-1)*dz(k-1))/(dz(k)+dz(k-1))!计算层界面电导率
enddo
do j=20,2000,20                     !j为厚度
    f=(500**2)*100/j                  !趋附深度求取频率
    E=(1d0,1d0)                     !迭代初值
    E(-n-1)=(1.0E-6,0d0)            !地上无穷远处电场1E-6v/m
    E(n+1)=(0d0,0d0)                !地下无穷远处电场为0
    do                              !SOR迭代法
      E1=E
      s=0
      do k=-n,n
            a=2d0*(E(k+1)/(dz(k)*dz(k)+dz(k-1)*dz(k))+E(k-1)/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1)))
            b=2d0/(dz(k)*dz(k)+dz(k-1)*dz(k))+2d0/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1))+i*f*miu*o1(k)
            E(k)=(1-w)*E(k)+w*a/b   
            s=s+abs(E(k)-E1(k))
      enddo
      if(s<1E-20)exit             !迭代结束条件
    enddo
    do k=-n-1,n
      H(k)=(E(k+1)-E(k))/(dz(k)*i*f*miu)      !计算各地层中点磁场
    enddo
    E14=0.75*E(0)+0.25d0*E(1)
    Hs=E14*o1(0)*dz(0)/2d0+H(0)               !地表处磁场
    ro=(abs(E(0)/Hs))**2/f/miu                  !卡尼亚视电阻率
    write(*,*)f,ro
    write(123,*)f,ro
enddo
pause
endprogram

aliouying 发表于 2014-4-21 10:21:50

本帖最后由 aliouying 于 2014-4-21 10:23 编辑

为什么要用差分呢？ sor迭代收敛很慢，一般需要迭代1000次以上，甚至几千次，你的程序有些问题，大致看了下，改了改program main
implicit none
integer,parameter::n=1000                     !地上n层 地下n层
complex*16,parameter::i=(0d0,1d0)               
real*8,parameter::PI=3.1415926535897932384626d0
real*8,parameter::miu=4.0d0*PI*1.0d-7            
real*8,parameter::w=0.8                        !SOR迭代因子
integer::j,k
real*8::o(-n-1:n),o1(-n:n),dz(-n-1:n+1),s,f,E14,ro,fa
!o地层电导率，o1分界面电导率，dz层厚度，f频率，E14四分之一厚度处电场，ro视电阻率，fa相位
complex*16::E(-n-1:n+1),H(-n-1:n),E1(-n-1:n+1),Hs,a,b
!E电场，H磁场，Hs地表磁场
open(123,file='a.dat')
dz=1d0                  !单层厚度1m
dz(-n-1)=1.0E6            !最上层厚度1000km，视作无限远
dz(n)=1.0E6               !最底层厚度1000km，视作无限远
forall (k=-n-1:-1:1) o(k)=1E-10      !空气中电导率1E-10，视作无限大
forall(k=0:300:1) o(k)=0.1d0         !0到300米，电导率0.1
forall(k=301:600:1) o(k)=0.005d0   !301到600米，电导率0.005
forall(k=601:n:1) o(k)=0.1d0         !601到无穷远，电导率0.1
do k=-n,n
    o1(k)=(o(k)*dz(k)+o(k-1)*dz(k-1))/(dz(k)+dz(k-1))!计算层界面电导率
enddo
do j=20,2000,20                     !j为厚度
    f=2*PI*(500**2)*100/j         !趋附深度求取频率
    E=(1d0,1d0)                     !迭代初值
    !E(-n-1)=(1.0E-6,0d0)         !地上无穷远处电场1E-6v/m
    E(n+1)=(0d0,0d0)                !地下无穷远处电场为0
    do                              !SOR迭代法
      E1=E
      s=0
      do k=-n,n
            a=2d0*(E(k+1)/(dz(k)*dz(k)+dz(k-1)*dz(k))+E(k-1)/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1)))
            b=2d0/(dz(k)*dz(k)+dz(k-1)*dz(k))+2d0/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1))+i*f*miu*o1(k)
            E(k)=(1-w)*E(k)+w*a/b   
            s=s+abs(E(k)-E1(k))**2
      enddo
      s = dsqrt( s/n/2 )
      if(s<1E-10)exit             !迭代结束条件
    enddo
    do k=-n-1,n
      H(k)=(E(k+1)-E(k))/(dz(k)*i*f*miu)      !计算各地层中点磁场
    enddo
    E14=0.75*E(0)+0.25d0*E(1)
    Hs=E14*o1(0)*dz(0)/2d0+H(0)               !地表处磁场
    ro=(abs(E(0)/Hs))**2/f/miu                  !卡尼亚视电阻率
    write(*,*)f,ro
    write(123,*)f,ro
enddo
pause
endprogram



geophylika 发表于 2014-5-20 00:04:00

做物探的不少，其实f 的代码现成的不少`
页: [1]
查看完整版本: MT一维正演，用SOR法不收敛