!----------------------------------------------------------------------- ! Laplace's equation: Gauss-Seidel method 伊理・藤野(1985) p.131 ! ! coded by Shigeo Yoden on Sept. 13, 2000 !----------------------------------------------------------------------- program laplace_gs use dcl ! DCL Fortran90 (dcl-5.1/dcl-f90) integer, parameter :: imax=20 real, parameter :: pi=3.141593, dx=pi/imax real, dimension(0:imax,0:imax) :: vold, vnew, van character :: ct*30, cl*30, csgi*1 !-- 境界条件 & 初期推定値 ---- do i=0, imax vold(i,imax) = sin(i*dx) end do do j=0, imax-1 vold(:,j) = j*dx*vold(:,imax) end do vold(imax,:) = 0 vnew = vold !!! write(*,*) 'epsilon=?' !!! read (*,*) epsilon epsilon = 0.0001 test = 2*epsilon lmn = 0 ! execute G-S iteration and calculate TEST where ! TEST is the maximum component of (vnew - vold) do while(test > epsilon) ! convergence check test = 0. do j=1, imax-1 do i=1, imax-1 vnew(i,j) = (vold(i-1,j) + vold(i,j-1) & + vold(i+1,j) + vold(i,j+1)) /4 diff = abs(vnew(i,j)-vold(i,j)) if(diff > test) test = diff end do end do vold = vnew lmn = lmn + 1 end do write(*,*) 'iteration # = ', lmn, ' DIFmax = ', test !-- 解析解 ---- do j=0, imax do i=0, imax van(i,j) = sinh(j*dx)*sin(i*dx) /sinh(pi) end do end do van(imax,:) = 0 !-- グラフ化 ---- call DclOpenGraphics() call DclNewFrame call DclSetWindow (0.0, 1.0, 0.0, 1.0) call DclSetViewPort(0.2, 0.8, 0.2, 0.8) call DclSetTransFunction call DclSetTitle('x', 'y', csgi(194)//' '//csgi(167), & csgi(194)//' '//csgi(167)) call DclDrawScaledAxis call DclSetParm('ENABLE_CONTOUR_MESSAGE', .false.) call DclSetParm('ENABLE_NEGATIVE_CONTOUR', .false.) call DclDrawContour(van) ! 解析解(実線) call DclSetParm('POSITIVE_CONTOUR_TYPE', 3) call DclDrawContour(vnew) ! 数値解(点線) ct = 'iteration # =?????' write(ct(14:18),'(i5)') lmn call DclDrawTitle('t', ct, position=-1.) cl = '=??.??E???s' write(cl(2:10),'(e9.2)') epsilon ct = csgi(156) // cl call DclDrawTitle('t', ct, position=-0.8) call DclDrawTitle('t', ' Gauss-Seidel method', position=-1.) call DclDrawTitle('t', 'Laplace'//csgi(178)//'s Equation:', & position=-1.) call DclCloseGraphics end program laplace_gs !----------------------------------------------------------------------- ! Laplace's equation: Successive overrelaxation(SOR) method ! Johnson and Riess, 1982: "Numerical Analysis(2nd edition)", ! Addison-Wesley Publishing Company, p.472 ! ! coded by Shigeo Yoden on Sept. 13, 2000 !----------------------------------------------------------------------- program laplace_sor use dcl ! DCL Fortran90 (dcl-5.1/dcl-f90) integer, parameter :: imax=20 real, parameter :: pi=3.141593, dx=pi/imax real, dimension(0:imax,0:imax) :: vold, vnew, van character :: ct*30, cl*30, csgi*1 !-- 境界条件 & 初期推定値 ---- do i=0, imax vold(i,imax) = sin(i*dx) end do do j=0, imax-1 vold(:,j) = j*dx*vold(:,imax) end do vold(imax,:) = 0 vnew = vold !!! write(*,*) 'epsilon=?' !!! read (*,*) epsilon epsilon = 0.001 test = 2*epsilon lmn = 0 !-- SOR ---- omega = 2/(1 + sin(dx)) ! execute SOR iteration and calculate TEST where ! TEST is the maximum component of (vnew - vold) do while(test > epsilon) ! convergence check test = 0. do j=1, imax-1 do i=1, imax-1 vnew(i,j) = vold(i,j) + omega*(vnew(i-1,j) + vnew(i,j-1) & + vold(i+1,j) + vold(i,j+1) -4*vold(i,j))/4 diff = abs(vnew(i,j)-vold(i,j)) if(diff > test) test = diff end do end do vold = vnew lmn = lmn + 1 end do write(*,*) 'iteration # = ', lmn, ' DIFmax = ', test !-- 解析解 ---- do j=0, imax do i=0, imax van(i,j) = sinh(j*dx)*sin(i*dx) /sinh(pi) end do end do van(imax,:) = 0 !-- グラフ化 ---- call DclOpenGraphics() call DclNewFrame call DclSetWindow (0.0, 1.0, 0.0, 1.0) call DclSetViewPort(0.2, 0.8, 0.2, 0.8) call DclSetTransFunction call DclSetTitle('x', 'y', csgi(194)//' '//csgi(167), & csgi(194)//' '//csgi(167)) call DclDrawScaledAxis call DclSetParm('ENABLE_CONTOUR_MESSAGE', .false.) call DclSetParm('ENABLE_NEGATIVE_CONTOUR', .false.) call DclDrawContour(van) ! 解析解(実線) call DclSetParm('POSITIVE_CONTOUR_TYPE', 3) call DclDrawContour(vnew) ! 数値解(点線) ct = 'iteration # =?????' write(ct(14:18),'(i5)') lmn call DclDrawTitle('t', ct, position=-1.) cl = '=??.??E???s' write(cl(2:10),'(e9.2)') epsilon ct = csgi(156) // cl call DclDrawTitle('t', ct, position=-0.8) call DclDrawTitle('t', ' SOR method', position=-1.) call DclDrawTitle('t', 'Laplace'//csgi(178)//'s Equation:', & position=-1.) call DclCloseGraphics end program laplace_sor !----------------------------------------------------------------------- ! Poisson's equation: Gauss-Seidel method 伊理・藤野(1985) p.131 ! ! coded by Shigeo Yoden on Sept. 13, 2000 !----------------------------------------------------------------------- program poisson_gs use dcl ! DCL Fortran90 (dcl-5.1/dcl-f90) integer, parameter :: imax=60 real, parameter :: pi=3.141593, dx=pi/imax real, dimension(0:imax,0:imax) :: vold, vnew, rho character :: ct*30, cl*30, csgi*1 !-- 境界条件 & 初期推定値 ---- vold = 0 vnew = 0 rho = 0 rho( imax/3, imax/3) = 1 rho(2*imax/3,2*imax/3) = - 2 !!! rho(imax/2,imax/2) = - 1 !!! write(*,*) 'epsilon=?' !!! read (*,*) epsilon epsilon = 0.000001 test = 2*epsilon lmn = 0 ! execute G-S iteration and calculate TEST where ! TEST is the maximum component of (vnew - vold) do while(test > epsilon) ! convergence check test = 0. do j=1, imax-1 do i=1, imax-1 vnew(i,j) = ( vnew(i-1,j) + vnew(i,j-1) & + vold(i+1,j) + vold(i,j+1) & - rho(i,j)*dx**2 ) /4 diff = abs(vnew(i,j)-vold(i,j)) if(diff > test) test = diff end do end do vold = vnew lmn = lmn + 1 end do write(*,*) 'iteration # = ', lmn, ' DIFmax = ', test !-- グラフ化 ---- call DclOpenGraphics() call DclNewFrame call DclSetWindow (0.0, 1.0, 0.0, 1.0) call DclSetViewPort(0.2, 0.8, 0.2, 0.8) call DclSetTransFunction call DclSetTitle('x', 'y', csgi(194)//' '//csgi(167), & csgi(194)//' '//csgi(167)) call DclDrawScaledAxis call DclDrawContour(vnew) ct = 'iteration # =?????' write(ct(14:18),'(i5)') lmn call DclDrawTitle('t', ct, position=-1.) cl = '=??.??E???s' write(cl(2:10),'(e9.2)') epsilon ct = csgi(156) // cl call DclDrawTitle('t', ct, position=-0.8) call DclDrawTitle('t', ' Gauss-Seidel method', position=-1.) call DclDrawTitle('t', 'Poisson'//csgi(178)//'s Equation:', position=-1.) call DclCloseGraphics end program poisson_gs