Fortran Coder

查看: 7948|回复: 9

[微分方程] 偏微分方程离散计算不收敛求助

[复制链接]

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
发表于 2014-5-13 16:12:15 | 显示全部楼层 |阅读模式
本帖最后由 岸边的鱼 于 2014-5-14 10:49 编辑

00000

00000

11111

11111

2222

2222

3333333

3333333

444444444

444444444

555555555555

555555555555

6666666666

6666666666

希望大家能给与指点,谢谢
我想使用Fortran语言编程求解如下一个偏微分方程组:

总结下我要求解的方程组即是上面的(4)(5)(6)式,我的FORTRAN程序是:

[Fortran] 纯文本查看 复制代码
module  BLSM
    implicit  none
integer,      parameter::     N=1000
real(kind=8), parameter::     dx=1000.d0
real(kind=8), parameter::     dy=1000.d0
real(kind=8), parameter::     dis=200000.d0
integer(kind=8), parameter::  M=-int(dis/dx)
real(kind=8), parameter::     tao=0.00036d0
real(kind=8), parameter::     v_F=5.556d0
real(kind=8), parameter::     h=10.d0
real(kind=8), parameter::     detp=00.d0/(1030.d0*9.8d0)
real(kind=8), parameter::     pi=3.1415926535898d0
real(kind=8), parameter::     f=0.0000773336d0
real(kind=8), parameter::     rho=1030.d0
real(kind=8), parameter::     g=9.8d0
real(kind=8), parameter::     ro=100000.d0
real(kind=8), parameter::     alpha=1E-10
real(kind=8), parameter::     lamda=1E-10
real(kind=8), parameter::     EPS=1E-28
endmodule BLSM
module YDHS
  use BLSM
  implicit none
contains
FunctionB_u(v,zeta,i,j,cishu)
implicitnone
real(kind=8)::B_u,P_Dx,P_Dx_ZH,P_Dx_DD,ii,jj,P_Dxx
real(kind=8)::v(M:N+M,-N-1:N+1)
real(kind=8)::zeta(M:N+M,-N:N),ii_2,jj_2
integer::i,j,cishu
     ii=i*dx
     jj=j*dy
     ii_2=ii*ii
     jj_2=jj*jj
B_u=-dy*f/(2.d0*V_F)*(v(i,j+1)+v(i-1,j+1)+v(i-1,j)+v(i,j))&
        +2.d0*g*dy/(V_F*dx)*(zeta(i,j)-zeta(i-1,j))&
        +(4.d0*dy*g*alpha*detp/V_F)*ii&
              *dexp(-alpha*ii_2-lamda*jj_2)&
        +(2.d0*dy*tao*ro**4.d0/(rho*h*V_F))/((ro*ro+ii_2+jj_2)**2.d0)&
              *(ii/2.d0+(dsqrt(3.d0)/2.d0)*jj)
return
endFunctionB_u

FunctionB_v(u,zeta,i,j,cishu)
implicitnone
real(kind=8)::B_v,P_Dy_ZH,ii,jj,ii_2,jj_2
real(kind=8)::u(M:N+M+1,-N:N)
real(kind=8)::zeta(M:N+M,-N:N)
integer::i,j,cishu,cishu_S_P_Dy
    ii=i*dx
    jj=j*dy
    ii_2=ii*ii
    jj_2=jj*jj
B_v=f*dy/(2.d0*V_F)*(u(i,j)+u(i,j-1)+u(i+1,j-1)+u(i+1,j))&
     +(2.d0*g/V_F)*(zeta(i,j)-zeta(i,j-1))&
+4.d0*g*lamda*detp*(dy/V_F)*jj*dexp(-alpha*ii_2-lamda*jj_2)&
-(2.d0*tao*ro**4.d0*dy/(rho*h*V_F))/((ro*ro+ii_2+jj_2)**2.0)&
                     *(dsqrt(3.d0)/2.d0*ii-jj/2.d0)
return
endFunctionB_v

Function  B_zeta(u,v,i,j)
implicit  none
real(kind=8)::B_zeta
real(kind=8)::u(M:N+M+1,-N:N)
real(kind=8)::v(M:N+M,-N-1:N+1)
integer::i,j
    B_zeta=(2.d0*h*dy/(V_F*dx))*(u(i+1,j)-u(i,j))&
       +(2.d0*h/V_F)*(v(i,j+1)-v(i,j))
return
end  Function  B_zeta
!______________
end  module YDHS

program  main
    use  BLSM
    use  YDHS
implicit  none
  real(kind=8)::ii,jj,P_Dx
real(kind=8)::u0(M:N+M+1,-N:N),v0(M:N+M,-N-1:N+1)
,zeta0(M:N+M,-N:N)
real(kind=8)::u1(M:N+M+1,-N:N),v1(M:N+M,-N-1:N+1)
,zeta1(M:N+M,-N:N)
  integer::i,j,cishu,k
  real(kind=8)::dx2
  integer::cishu_u,cishu_v,cishu_zeta
  character*80::filename
  open (30,file='P_Dy_ZH.dat')
  open (20,file='P_Dx_ZH.dat')
!初值全部设定为0
    cishu=0
    u0=0.d0
    v0=0.d0
    zeta0=0.d0
!  先赋初值再使用
       u1=0.d0
       v1=0.d0
       zeta1=0.d0
!-----------------------
      u1(M,:)  =0.d0      !u的岸边界值
      u1(M+1,:) =0.d0      !u的第一次迭代边界值
      v1(:,-N) =0.d0      !v的水边界值
      v1(:,-N+1)=0.d0      !v第一次迭代的边界值
      zeta1(:,:)=0.d0
100 write(*,*)  "迭代次数",cishu
        cishu=cishu+1
        cishu_u=cishu+10000
       write(filename,FMT='(I6.6,A4)')cishu_u,'.dat'
            filename=trim(filename)
open(100,file=filename,status='unknown')
       doi=M+1,N+M
          ii=i*dx
           if(i==300.or.i==600)then
               write(*,*) i,'I''m_u'
           end if
                u1(i,-N)=u1(i,-N+1)   !输运条件导数为0
              doj=-N+1,N-1
                 jj=j*dy
                u1(i,j+1)=u1(i,j-1)+B_u(v0,zeta0,i,j,cishu)
                    write(100,*)ii,jj, u1(i,j+1)
               end do
         end  do
close (100)
!----------------------
        cishu_v=cishu+20000
   write(filename,FMT='(I6.6,A4)')cishu_v,'.dat'
         filename=trim(filename)
       open(101,file=filename,status='unknown')
   do  i=M,N+M
       ii=i*dx
       if(i==300.or.i==600)then
         write(*,*) i,'I''m_v'
       end if
           v1(:,-N-1)=v1(:,-N)   !提输运条件
       do  j=-N+1,N
           jj=j*dy
          v1(i,j+1)=v1(i,j-1)+B_v(u0,zeta0,i,j,cishu)
           write(101,*)   ii,jj, v1(i,j+1)
end do
   end do
close (101)
!------------------------------------------------
cishu_zeta=cishu+30000
write(filename,FMT='(I6.6,A4)')cishu_zeta,'.dat'
      filename=trim(filename)
open(1001,file=filename,status='unknown')
      do  i=M,N+M
          ii=i*dx
          if(i==300.or.i==600)then
             write(*,*) i,'Im_zeta'
          end if
        doj=-N+1,N-1
           jj=j*dy
          zeta1(i,j+1)=zeta1(i,j-1)+B_zeta(u1,v1,i,j)
           write(1001,*)  ii,jj, zeta1(i,j+1)   !030001.dat
         end do
      end do
close (1001)
!________判断前后两次差别[/align]

[align=left]  do i=M,N+M
    do j=-N,N
dx2=dx2+(u1(i,j)-u0(i,j))**2.d0+(v1(i,j)-v0(i,j))**2.d0
+(zeta1(i,j)-zeta0(i,j))**2.d0
    end do
  end do
       dx2=dsqrt(dx2)
write(*,*) "前后差距",dx2
    if(dx2>=eps)then
       u0=u1
       v0=v1
       zeta0=zeta1
       goto  100
    end if
stop
endprogrammain

计算显示计算不收敛。。求助各位大神帮忙看下是程序的问题还是本身这种方法就是不收敛的?

科研穷三代,读博毁一生

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
 楼主| 发表于 2014-5-13 16:14:02 | 显示全部楼层
因为长度受到限制,所以就没有吧计算结果的图发上来,大家可以运行一下,迭代下去会逐步变大的
科研穷三代,读博毁一生

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
 楼主| 发表于 2014-5-13 16:14:40 | 显示全部楼层
感谢各位,帮忙看看的大神,这个破问题纠结了2个月了,郁闷死了
科研穷三代,读博毁一生

953

帖子

0

主题

0

精华

大师

F 币
180 元
贡献
73 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
发表于 2014-5-13 22:37:27 | 显示全部楼层
太长了,实在没精力看。
我觉得应该是楼主的代码写得不好,而不是算法问题。

你确定 dx2 就只算最后一个?而不是累加?

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
 楼主| 发表于 2014-5-14 10:10:38 | 显示全部楼层
vvt 发表于 2014-5-13 22:37
太长了,实在没精力看。
我觉得应该是楼主的代码写得不好,而不是算法问题。

我的dx2就是累加啊?不是算最后一个,dx2是u1,v1,zeta1 三个数组对应的数字与前面一次计算的u0,v0,zeta0的差值啊。还是谢谢你的关注。
科研穷三代,读博毁一生

953

帖子

0

主题

0

精华

大师

F 币
180 元
贡献
73 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
发表于 2014-5-14 10:17:29 | 显示全部楼层
[Fortran] 纯文本查看 复制代码
do i=M,N+M
do j=-N,N
  dx2=(u1(i,j)-u0(i,j))**2.d0+(v1(i,j)-v0(i,j))**2.d0+(zeta1(i,j)-zeta0(i,j))**2.d0
  end do
end do


这就是只算了最后一个。不是累加。累加应该是 dx2 = dx2 + ...

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
 楼主| 发表于 2014-5-14 10:20:29 | 显示全部楼层
我的神啊,怎么犯了这么个低级错误啊,我要头撞南墙了,无地自容啊,跪谢大神
科研穷三代,读博毁一生

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
 楼主| 发表于 2014-5-14 10:44:03 | 显示全部楼层
本帖最后由 岸边的鱼 于 2014-5-14 10:51 编辑

感谢楼上的大神,指出这个错误,不幸的是我修改后还是不收敛,原帖中的这个错误我已经在原帖中修改了再次感谢。 QQ图片20140514104345.jpg
科研穷三代,读博毁一生

136

帖子

3

主题

0

精华

版主

F 币
1964 元
贡献
1677 点

帅哥勋章管理勋章爱心勋章新人勋章热心勋章元老勋章

发表于 2014-5-16 17:02:15 | 显示全部楼层
算法看得不是很明白,最近毕业的事情都快疯了~  没太多时间帮忙看

66

帖子

5

主题

2

精华

版主

院士级水师

F 币
481 元
贡献
273 点

管理勋章帅哥勋章爱心勋章规矩勋章

QQ
 楼主| 发表于 2014-5-16 19:47:18 | 显示全部楼层
aliouying 发表于 2014-5-16 17:02
算法看得不是很明白,最近毕业的事情都快疯了~  没太多时间帮忙看

好的,谢谢热心关注,祝你答辩顺利
科研穷三代,读博毁一生
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-2-25 21:02

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表