! 拡散方程式の陰解法(backword 法) ! \DP{u}{t} = \kappa \DP[2]{u}{x} ! \kappa=0.01 ! program back implicit none real(8), parameter :: pi = 4.0d0*atan(1.0d0) real(8), parameter :: kappa = 0.01d0 real(8), parameter :: xmin = 0.0d0 real(8), parameter :: xmax = 1.0d0 real(8), parameter :: tmin = 0.0d0 real(8), parameter :: tmax = 10.0d0 real(8), allocatable :: U_a(:) !U^{n+1} real(8), allocatable :: U_b(:) !U^{n} real(8), allocatable :: X_cal(:) !xの値 real(8), allocatable :: L_mat(:,:) real(8), allocatable :: U_mat(:,:) real(8), allocatable :: A_mat(:,:) real, allocatable :: U_sing(:) real, allocatable :: X_sing(:) real, allocatable :: genmitu_sing(:) real(8) :: dx real(8) :: dt real(8) :: f_i real(8) :: LU integer :: x_step integer :: t_step integer :: x_lupe integer :: t_lupe integer :: i_lupe integer :: k_lupe integer :: j_lupe integer :: IWS integer :: scheme character(10) chara write(*,*) "Enter in time steps. & & if you want to enter more 10000000000,& & enter in 'a'. and rekey. default is 800."! 500. read(*,'(a)') chara !時間分割数の仮入力 call empty(chara, 200, t_step) !仮文字列が空白の場合デフォルト値200を入れる write(*,*) "You chose ", t_step write(*,*) "Enter in X steps. & & if you want to enter more 10000000000,& & enter in 'a'. and rekey. default is 10."!20 " read(*,'(a)') chara !空間分割数の仮文字列入力 call empty(chara, 50, x_step) ! call empty(chara, 3, x_step) !仮文字列が空白の場合デフォルト値10を入れる write(*,*) "You chose", x_step open(10,file='test.dat', status='replace') allocate(U_a(0:x_step)) allocate(U_b(0:x_step)) allocate(X_cal(0:x_step)) allocate(U_sing(0:x_step)) allocate(X_sing(0:x_step)) allocate(A_mat(1:x_step-1,1:x_step-1)) allocate(L_mat(1:x_step-1,1:x_step-1)) allocate(U_mat(1:x_step-1,1:x_step-1)) allocate(genmitu_sing(0:x_step)) dx = (abs(xmax) - abs(xmin))/dble(x_step) dt = (abs(tmax) - abs(tmin))/dble(t_step) !------------------------------------------ call swpset('LDUMP',.true.) call swpset('LWAIT',.false.) call swpset('lwait1',.false.) call sgpwsn read(*,'(a)') chara call empty(chara,1,IWS) call gropn(IWS) !-----初期条件---------------------------------- do x_lupe = 0, x_step X_cal(x_lupe) = xmin + dx*dble(x_lupe) U_b(x_lupe) = sin(pi*X_cal(x_lupe)) end do X_sing=real(X_cal) U_sing=real(U_b) call grfrm call grswnd(0.0,1.0,0.0,1.0) call uspfit call grstrf call usdaxs call uulin(x_step+1, X_sing(0:x_step), U_sing(0:x_step)) A_mat(:,:)=0.0d0 U_mat(:,:)=0.0d0 L_mat(:,:)=0.0d0 LU = 0.0d0 U_b(0) = U_b(x_step) U_a = U_b do x_lupe = 1, x_step-1 A_mat(x_lupe,x_lupe) = 1.0d0 + 2.0d0*(kappa*dt/(dx**2)) A_mat(x_lupe,x_lupe+1) = - (kappa*dt/(dx**2)) A_mat(x_lupe+1,x_lupe) = - (kappa*dt/(dx**2)) end do ! A_mat(x_step-1,x_step-1) = 1.0d0 + 2.0d0*(kappa*dt/(dx**2)) ! A_mat(1,1) = 1.0d0 ! A_mat(1,2) = 2.0d0 ! A_mat(1,3) = 3.0d0 ! A_mat(1,4) = 4.0d0 ! A_mat(2,1) = 2.0d0 ! A_mat(2,2) = 3.0d0 ! A_mat(2,3) = 4.0d0 ! A_mat(2,4) = 1.0d0 ! A_mat(3,1) = 3.0d0 ! A_mat(3,2) = 4.0d0 ! A_mat(3,3) = 1.0d0 ! A_mat(3,4) = 2.0d0 ! A_mat(4,1) = 4.0d0 ! A_mat(4,2) = 1.0d0 ! A_mat(4,3) = 2.0d0 ! A_mat(4,4) = 3.0d0 ! A_mat(1,1) = 1.0d0 ! A_mat(1,2) = 2.0d0 ! A_mat(1,3) = 3.0d0 ! A_mat(2,1) = 2.0d0 ! A_mat(2,2) = 1.0d0 ! A_mat(2,3) = 2.0d0 ! A_mat(3,1) = 3.0d0 ! A_mat(3,2) = 2.0d0 ! A_mat(3,3) = 1.0d0 !-----LU分解-------------------------------- ! 1行目の計算 ! do i_lupe = 1,x_step ! ! L_mat(i_lupe,1) = A_mat(i_lupe, 1) ! ! U_mat(1,i_lupe) = A_mat(1,i_lupe)/L_mat(1,1) ! end do ! ! !2行目以降 ! do k_lupe = 2, x_step - 1 ! ! L_mat(k_lupe, 2) = & ! & A_mat(k_lupe,2) - L_mat(k_lupe,1) * U_mat(1,2) ! ! do i_lupe = 3, k_lupe ! ! do j_lupe = 2, i_lupe-1 ! ! LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe) ! ! end do ! ! L_mat(k_lupe,i_lupe) = & ! & A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe) & ! & - LU ! ! LU = 0.0d0 ! ! end do ! ! ! do i_lupe = k_lupe + 1, x_step ! ! do j_lupe = 2, k_lupe -1 ! ! LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe) ! end do ! ! U_mat(k_lupe,i_lupe) = & ! & (A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe)& ! & - LU)/L_mat(k_lupe,k_lupe) ! ! end do ! ! end do ! ! LU = 0.0d0 ! !x_step行目の計算 ! do j_lupe = 2, x_step-1 ! ! LU = LU + L_mat(x_step,j_lupe)*U_mat(j_lupe,x_step) ! ! end do ! L_mat(x_step,x_step) = & ! & A_mat(x_step,x_step) - L_mat(x_step,1)*U_mat(1,x_step) & ! & - LU ! ! do x_lupe = 1, x_step ! do t_lupe = 1, x_step ! write(10,*) (L_mat(x_lupe,t_lupe)),U_mat(x_lupe,t_lupe) ! end do ! end do call LUdecom(x_step-1,A_mat,L_mat,U_mat) ! do x_lupe = 1, x_step ! do t_lupe = 1, x_step ! write(10,*) L_mat(x_lupe,t_lupe),U_mat(x_lupe,t_lupe) ! end do ! end do ! write(10,*) "(1,1)", L_mat(1,1)*U_mat(1,1), & ! & "(1,2)", L_mat(1,1)*U_mat(1,2), & ! & "(1,3)", L_mat(1,1)*U_mat(1,3) ! ! write(10,*) "(2,1)", L_mat(2,1)*U_mat(1,1), & ! & "(2,2)", L_mat(2,1)*U_mat(1,2) + L_mat(2,2), & ! & "(2,3)", L_mat(2,1)*U_mat(1,3) + L_mat(2,2)*U_mat(2,3) ! ! write(10,*) "(3,1)", L_mat(3,1)*U_mat(1,1), & ! & "(3,2)", L_mat(3,1)*U_mat(1,2) + L_mat(3,2), & ! & "(3,3)", L_mat(3,1)*U_mat(1,3) + L_mat(3,2)*U_mat(2,3) +& ! & L_mat(3,3) ! write(10,*) "(1,1)", L_mat(1,1)*U_mat(1,1), & ! & "(1,2)", L_mat(1,1)*U_mat(1,2), & ! & "(1,3)", L_mat(1,1)*U_mat(1,3), & ! & "(1,4)", L_mat(1,1)*U_mat(1,4) ! write(10,*) "(2,1)", L_mat(2,1)*U_mat(1,1), & ! & "(2,2)", L_mat(2,1)*U_mat(1,2) + L_mat(2,2), & ! & "(2,3)", L_mat(2,1)*U_mat(1,3) + L_mat(2,2)*U_mat(2,3), & ! & "(2,4)", L_mat(2,1)*U_mat(1,4) + L_mat(2,2)*U_mat(2,4) ! write(10,*) "(3,1)", L_mat(3,1)*U_mat(1,1), & ! & "(3,2)", L_mat(3,1)*U_mat(1,2) + L_mat(3,2), & ! & "(3,3)", L_mat(3,1)*U_mat(1,3) + L_mat(3,2)*U_mat(2,3) +& ! & L_mat(3,3), & ! & "(3,4)", L_mat(3,1)*U_mat(1,4) + L_mat(3,2)*U_mat(2,4) +& ! & L_mat(3,3)*U_mat(3,4) ! write(10,*) "(4,1)", L_mat(4,1)*U_mat(1,1), & ! & "(4,2)", L_mat(4,1)*U_mat(1,2) + L_mat(4,2), & ! & "(4,3)", L_mat(4,1)*U_mat(1,3) + L_mat(4,2)*U_mat(2,3) +& ! & L_mat(4,3), & ! & "(4,4)", L_mat(4,1)*U_mat(1,4) + L_mat(4,2)*U_mat(2,4) +& ! & L_mat(4,3)*U_mat(3,4) + L_mat(4,4) ! do x_lupe = 1,3 ! ! write(10,*) "(",x_lupe,"1)", L_mat(x_lupe,1), & ! & "(",x_lupe,"2)", L_mat(x_lupe,2), & ! & "(",x_lupe,"3)", L_mat(x_lupe,3) ! ! end do ! write(10,*) "" ! do x_lupe = 1,3 ! ! write(10,*) "(",x_lupe,"1)", U_mat(x_lupe,1), & ! & "(",x_lupe,"2)", U_mat(x_lupe,2), & ! & "(",x_lupe,"3)", U_mat(x_lupe,3) ! ! end do do t_lupe = 1, t_step !xyの計算 U_a(1) = U_b(1)/L_mat(1,1) do x_lupe = 2, x_step-1 LU = 0.0d0 do j_lupe = 1, x_lupe - 1 LU = LU + L_mat(x_lupe,j_lupe) * U_a(j_lupe) end do U_a(x_lupe) = (U_b(x_lupe) - LU)/L_mat(x_lupe,x_lupe) end do ! U_b(1:x_step-1) = U_a(x_step-1:1:-1) ! do x_lupe = 2, x_step-1 ! ! LU = 0.0d0 ! ! do j_lupe = 1, x_lupe-1 ! ! LU = U_mat(x_lupe,j_lupe) * U_b(j_lupe) ! ! end do ! ! U_b(x_lupe) = U_a(x_lupe) - LU ! ! end do do x_lupe = x_step-1, 1, -1 LU = 0.0d0 do j_lupe = x_step-1, x_lupe+1, -1 LU = U_mat(x_lupe,j_lupe) * U_b(j_lupe) end do U_b(x_lupe) = U_a(x_lupe) - LU end do U_b(0) = 0.0d0 U_b(x_step) = 0.0d0 ! U_a(0) = U_b(0) ! U_a(x_step) = U_b(x_step) ! U_a(1:x_step-1) = U_b(x_step-1:1:-1) ! U_sing = real(U_a) U_sing = real(U_b) ! U_b = U_a do x_lupe = 0, x_step genmitu_sing(x_lupe) = & & sin(pi*X_sing(x_lupe))*exp(-(pi**2)*kappa*dt*t_lupe) end do call grfrm call grswnd(0.0,1.0,0.0,1.0) call usspnt(x_step+1,X_sing(0:x_step),U_sing(0:x_step)) call usspnt(x_step+1,X_sing(0:x_step),genmitu_sing(0:x_step)) call uspfit call grstrf call usdaxs call uuslnt(1) call uulin(x_step+1, X_sing(0:x_step), U_sing(0:x_step)) call uuslnt(2) call uulin(x_step+1, X_sing(0:x_step), genmitu_sing(0:x_step)) end do call grcls close(10) deallocate(U_a,U_b,X_sing,X_cal,U_sing,L_mat,A_mat,U_mat) stop end program back !-----LU分解------------------------------------- subroutine LUdecom(x_step, A_mat,L_mat,U_mat) implicit none integer k_lupe integer i_lupe integer j_lupe integer, intent(in) :: x_step real(8), intent(in) :: A_mat(1:x_step,1:x_step) real(8), intent(out) :: L_mat(1:x_step,1:x_step) real(8), intent(out) :: U_mat(1:x_step,1:x_step) real(8) :: LU LU = 0.0d0 ! do k_lupe = 1, x_step ! ! U_mat(k_lupe,k_lupe) = 1.0d0 ! ! end do do i_lupe = 1,x_step L_mat(i_lupe,1) = A_mat(i_lupe, 1) U_mat(1,i_lupe) = A_mat(1,i_lupe)/L_mat(1,1) end do !2行目以降 do k_lupe = 2, x_step L_mat(k_lupe, 2) = & & A_mat(k_lupe,2) - L_mat(k_lupe,1) * U_mat(1,2) ! U_mat(k_lupe, 3) = & ! & (A_mat(k_lupe, 3) - L_mat(k_lupe,1)*U_mat(1,3))/ l_mat( do i_lupe = 3, k_lupe do j_lupe = 2, i_lupe-1 LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe) end do L_mat(k_lupe,i_lupe) = & & A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe) & & - LU LU = 0.0d0 end do do i_lupe = k_lupe + 1, x_step do j_lupe = 2, k_lupe -1 LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe) end do U_mat(k_lupe,i_lupe) = & & (A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe)& & - LU)/L_mat(k_lupe,k_lupe) LU = 0.0d0 end do end do ! LU = 0.0d0 !x_step行目の計算 ! do j_lupe = 2, x_step -1 ! ! LU = LU + L_mat(x_step,j_lupe)*U_mat(j_lupe,x_step) ! ! end do ! ! L_mat(x_step,x_step) = & ! & A_mat(x_step,x_step) - L_mat(x_step,1)*U_mat(1,x_step) & ! & - LU ! return end subroutine LUdecom !-----入力判断サブルーチン------------------------- !入力診断のサブルーチン. !キーボード入力が空白のときはデフォルト値を入れる ! subroutine empty(chara,default_step,step) !------入力診断の宣言文---------------------- implicit none character(10) chara integer default_step integer step !------入力診断------------------------------ if(len_trim(chara) == 0)then step = default_step !空白の場合デフォルト値を代入 else if(chara == "a")then read(*,*) step !aの場合さらにもう一度代入 else read(chara,*) step !仮文字列から実際の変数への代入 end if end if return end subroutine empty