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

岸边的鱼

希望大家能给与指点，谢谢

我想使用Fortran语言编程求解如下一个偏微分方程组：

总结下我要求解的方程组即是上面的（4）（5）（6）式，我的FORTRAN程序是：

[Fortran] 纯文本查看 复制代码

查看代码

打印?

001	module  BLSM

002	implicit  none

003	integer,      parameter::     N=1000

004	real(kind=8), parameter::     dx=1000.d0

005	real(kind=8), parameter::     dy=1000.d0

006	real(kind=8), parameter::     dis=200000.d0

007	integer(kind=8), parameter::  M=-int(dis/dx)

008	real(kind=8), parameter::     tao=0.00036d0

009	real(kind=8), parameter::     v_F=5.556d0

010	real(kind=8), parameter::     h=10.d0

011	real(kind=8), parameter::     detp=00.d0/(1030.d0*9.8d0)

012	real(kind=8), parameter::     pi=3.1415926535898d0

013	real(kind=8), parameter::     f=0.0000773336d0

014	real(kind=8), parameter::     rho=1030.d0

015	real(kind=8), parameter::     g=9.8d0

016	real(kind=8), parameter::     ro=100000.d0

017	real(kind=8), parameter::     alpha=1E-10

018	real(kind=8), parameter::     lamda=1E-10

019	real(kind=8), parameter::     EPS=1E-28

020	endmodule BLSM

021	module YDHS

022

use BLSM

023	implicit none

024

contains

025	FunctionB_u(v,zeta,i,j,cishu)

026	implicitnone

027	real(kind=8)::B_u,P_Dx,P_Dx_ZH,P_Dx_DD,ii,jj,P_Dxx

028	real(kind=8)::v(M:N+M,-N-1:N+1)

029	real(kind=8)::zeta(M:N+M,-N:N),ii_2,jj_2

030	integer::i,j,cishu

031

ii=i*dx

032

jj=j*dy

033	ii_2=ii*ii

034	jj_2=jj*jj

035	B_u=-dy*f/(2.d0*V_F)*(v(i,j+1)+v(i-1,j+1)+v(i-1,j)+v(i,j))&

036	+2.d0*g*dy/(V_F*dx)*(zeta(i,j)-zeta(i-1,j))&

037	+(4.d0*dy*g*alpha*detp/V_F)*ii&

038	*dexp(-alpha*ii_2-lamda*jj_2)&

039	+(2.d0*dy*tao*ro**4.d0/(rho*h*V_F))/((ro*ro+ii_2+jj_2)**2.d0)&

040	*(ii/2.d0+(dsqrt(3.d0)/2.d0)*jj)

041

return

042	endFunctionB_u

043

044	FunctionB_v(u,zeta,i,j,cishu)

045	implicitnone

046	real(kind=8)::B_v,P_Dy_ZH,ii,jj,ii_2,jj_2

047	real(kind=8)::u(M:N+M+1,-N:N)

048	real(kind=8)::zeta(M:N+M,-N:N)

049	integer::i,j,cishu,cishu_S_P_Dy

050

ii=i*dx

051

jj=j*dy

052	ii_2=ii*ii

053	jj_2=jj*jj

054	B_v=f*dy/(2.d0*V_F)*(u(i,j)+u(i,j-1)+u(i+1,j-1)+u(i+1,j))&

055	+(2.d0*g/V_F)*(zeta(i,j)-zeta(i,j-1))&

056	+4.d0*g*lamda*detp*(dy/V_F)*jj*dexp(-alpha*ii_2-lamda*jj_2)&

057	-(2.d0*tao*ro**4.d0*dy/(rho*h*V_F))/((ro*ro+ii_2+jj_2)**2.0)&

058	*(dsqrt(3.d0)/2.d0*ii-jj/2.d0)

059

return

060	endFunctionB_v

061

062	Function  B_zeta(u,v,i,j)

063	implicit  none

064	real(kind=8)::B_zeta

065	real(kind=8)::u(M:N+M+1,-N:N)

066	real(kind=8)::v(M:N+M,-N-1:N+1)

067	integer::i,j

068	B_zeta=(2.d0*h*dy/(V_F*dx))*(u(i+1,j)-u(i,j))&

069	+(2.d0*h/V_F)*(v(i,j+1)-v(i,j))

070

return

071	end  Function  B_zeta

072	!______________

073	end  module YDHS

074

075	program  main

076

use  BLSM

077

use  YDHS

078	implicit  none

079	real(kind=8)::ii,jj,P_Dx

080	real(kind=8)::u0(M:N+M+1,-N:N),v0(M:N+M,-N-1:N+1)

081	,zeta0(M:N+M,-N:N)

082	real(kind=8)::u1(M:N+M+1,-N:N),v1(M:N+M,-N-1:N+1)

083	,zeta1(M:N+M,-N:N)

084	integer::i,j,cishu,k

085	real(kind=8)::dx2

086	integer::cishu_u,cishu_v,cishu_zeta

087	character*80::filename

088	open (30,file='P_Dy_ZH.dat')

089	open (20,file='P_Dx_ZH.dat')

090

!初值全部设定为0

091

cishu=0

092

u0=0.d0

093

v0=0.d0

094	zeta0=0.d0

095	!  先赋初值再使用

096

u1=0.d0

097

v1=0.d0

098	zeta1=0.d0

099	!-----------------------

100	u1(M,:)  =0.d0      !u的岸边界值

101	u1(M+1,:) =0.d0      !u的第一次迭代边界值

102	v1(:,-N) =0.d0      !v的水边界值

103	v1(:,-N+1)=0.d0      !v第一次迭代的边界值

104	zeta1(:,:)=0.d0

105	100 write(*,*)  "迭代次数",cishu

106	cishu=cishu+1

107	cishu_u=cishu+10000

108	write(filename,FMT='(I6.6,A4)')cishu_u,'.dat'

109	filename=trim(filename)

110	open(100,file=filename,status='unknown')

111	doi=M+1,N+M

112

ii=i*dx

113	if(i==300.or.i==600)then

114	write(*,*) i,'I''m_u'

115

end if

116	u1(i,-N)=u1(i,-N+1)   !输运条件导数为0

117	doj=-N+1,N-1

118

jj=j*dy

119	u1(i,j+1)=u1(i,j-1)+B_u(v0,zeta0,i,j,cishu)

120	write(100,*)ii,jj, u1(i,j+1)

121

end do

122

end  do

123	close (100)

124	!----------------------

125	cishu_v=cishu+20000

126	write(filename,FMT='(I6.6,A4)')cishu_v,'.dat'

127	filename=trim(filename)

128	open(101,file=filename,status='unknown')

129	do  i=M,N+M

130

ii=i*dx

131	if(i==300.or.i==600)then

132	write(*,*) i,'I''m_v'

133

end if

134	v1(:,-N-1)=v1(:,-N)   !提输运条件

135	do  j=-N+1,N

136

jj=j*dy

137	v1(i,j+1)=v1(i,j-1)+B_v(u0,zeta0,i,j,cishu)

138	write(101,*)   ii,jj, v1(i,j+1)

139

end do

140

end do

141	close (101)

142	!------------------------------------------------

143	cishu_zeta=cishu+30000

144	write(filename,FMT='(I6.6,A4)')cishu_zeta,'.dat'

145	filename=trim(filename)

146	open(1001,file=filename,status='unknown')

147	do  i=M,N+M

148

ii=i*dx

149	if(i==300.or.i==600)then

150	write(*,*) i,'Im_zeta'

151

end if

152	doj=-N+1,N-1

153

jj=j*dy

154	zeta1(i,j+1)=zeta1(i,j-1)+B_zeta(u1,v1,i,j)

155	write(1001,*)  ii,jj, zeta1(i,j+1)   !030001.dat

156

end do

157

end do

158	close (1001)

159	!________判断前后两次差别[/align]

160

161	[align=left]  do i=M,N+M

162

do j=-N,N

163	dx2=dx2+(u1(i,j)-u0(i,j))**2.d0+(v1(i,j)-v0(i,j))**2.d0

164	+(zeta1(i,j)-zeta0(i,j))**2.d0

165

end do

166

end do

167	dx2=dsqrt(dx2)

168	write(*,*) "前后差距",dx2

169	if(dx2>=eps)then

170

u0=u1

171

v0=v1

172	zeta0=zeta1

173

goto  100

174

end if

175

stop

176	endprogrammain

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

岸边的鱼

因为长度受到限制，所以就没有吧计算结果的图发上来，大家可以运行一下，迭代下去会逐步变大的

岸边的鱼

感谢各位，帮忙看看的大神，这个破问题纠结了2个月了，郁闷死了

vvt

太长了，实在没精力看。
我觉得应该是楼主的代码写得不好，而不是算法问题。

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

岸边的鱼

vvt 发表于 2014-5-13 22:37
太长了，实在没精力看。
我觉得应该是楼主的代码写得不好，而不是算法问题。

我的dx2就是累加啊？不是算最后一个，dx2是u1，v1，zeta1 三个数组对应的数字与前面一次计算的u0，v0，zeta0的差值啊。还是谢谢你的关注。

vvt

[Fortran] 纯文本查看 复制代码

查看代码

打印?

1	do i=M,N+M

2

do j=-N,N

3	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

4

end do

5

end do

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

岸边的鱼

我的神啊，怎么犯了这么个低级错误啊，我要头撞南墙了，无地自容啊，跪谢大神

岸边的鱼

感谢楼上的大神，指出这个错误，不幸的是我修改后还是不收敛,原帖中的这个错误我已经在原帖中修改了再次感谢。

aliouying

算法看得不是很明白，最近毕业的事情都快疯了~  没太多时间帮忙看

岸边的鱼

aliouying 发表于 2014-5-16 17:02
算法看得不是很明白，最近毕业的事情都快疯了~  没太多时间帮忙看

好的，谢谢热心关注，祝你答辩顺利