请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 1927|回复: 0

[求助] 关于张弛法编程的问题

[复制链接]

新浪微博达人勋

发表于 2015-8-10 22:51:15 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
下面给的代码是我写的顺序张弛求解椭圆方程的代码,当超张弛系数(代码里c)为1的时候可以正常出结果(程序会运行10多秒),但是系数一旦比1大一点点,立马程序就会结束运行,而且出了结果也不对,我试过比1小的系数,也是可以出结果的,但是时间比系数是1的长,但是书上说这个系数在1~2之间啊,有没有大神帮忙解决下!!!

------------------------------------------------------------------------------------------------------------------------------------------------
program main
implicit none
!声明变量与常数
integer::irec,i,j,flag,a,b
integer,parameter::m=101,n=1001
real,parameter::pi=3.14159,d=pi/50.0,w=0.0001,c=1
real::g(m,n),r(m,n)
real::start,finish
call cpu_time(start)

!给上下边界赋值
do i=1,m
   g(i,1)=-sin((i-1)*d)
   g(i,n)=-sin((i-1)*d)*cos((n-1)*d)
end do

!给定侧边界值
do j=2,n
    g(1,j)=0
    g(101,j)=0
    end do

!给估计值
do j=2,n-1
    do i=1,m-1
        g(i,j)=-sin((i-1)*d)
    end do
end do

flag=1

!求helmholtz方程的近似解

        a=2
        b=2

        !利用同时张弛法进行迭代逼近
        do while(flag==1)
        
        !同时算第m次所有点的余差
        do j=2,n-1
            do i=2,m-1
                r(i,j)=g(i-1,j)+g(i+1,j)+g(i,j-1)+g(i,j+1)-(4+2*d*d)*g(i,j)-4*d*d*sin((i-1)*d)*cos((j-1)*d)
            end do
        end do
            !对所有点的余差r进行判断
             flag=0
            loop1:do j=2,n-1
                do i=2,m-1
                    if(abs(r(i,j)/(4+2*d*d))>w) then                        
                        flag=1
                        a=i
                        b=j      
                        exit loop1
                    end if
                end do               
            end do loop1

          !求算m+1次的估计值
            do j=b,n-1
                do i=a,m-1                  
                 g(i,j)=g(i,j)+c*r(i,j)/(4+2*d*d)  
                end do
            end do
            
        end do
          call cpu_time(finish)
write(*,*) finish-start
   !写.grd文件
   open(unit=11,file="G:\g(x,y)1.grd",form="unformatted",access="direct",status="unknown",recl=1)
        irec=1
   do j=1,n
            do i=1,m
                                write(11,rec=irec) g(i,j)
                        irec=irec+1
                end do
    end do
    close(11)

end program
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

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