#2
by 太郎 » 5年前
コードを貼るの忘れてました
コード:
program main
implicit none
double precision fai(10000),t,dt,x,dx,xmax
dx=0.1d0
xmax=2.0d0
dt=0.001d0
t=0.05d0
call youkaihou(fai,t,dt,x,dx,xmax)
call inkaihou(fai,t,dt,x,dx,xmax)
!call cranknicholson(fai,t,dt,x,dx,xmax)
end program main
subroutine youkaihou(fai,t,dt,x,dx,xmax)!t秒後の温度を求める。
implicit none
double precision fai0(10000),fai(10000),dum(10000),t,dt,x,dx,xmax,pi,r
integer n,m,i,j
pi=3.14159265359d0
x=0.0d0
n=dble(xmax/dx)
do i=1,n
x=x+dx
fai0(i)=sin(pi*x/2.0d0)
enddo
m=dble(t/dt)
r=dt/(dx**2.0d0)
do i=1,m
do j=1,n
if(i==1)then
if(j==1) then
dum(j)=(1.0d0-2.0d0*r)*fai0(j)+r*fai0(j+1)
else if(j==n) then
dum(j)=r*fai0(n-1)+(1.0d0-2.0d0*r)*fai0(n)
else!1<j<n-1のとき
dum(j)=r*fai0(j-1)+(1.0d0-2.0d0*r)*fai0(j)+r*fai0(j)
endif
else !i>2のとき
if(j==1) then
dum(j)=(1.0d0-2.0d0*r)*fai(j)+r*fai(j+1)
else if(j==n) then
dum(j)=r*fai(n-1)+(1.0d0-2.0d0*r)*fai(n)
else!1<j<n-1のとき
dum(j)=r*fai(j-1)+(1.0d0-2.0d0*r)*fai(j)+r*fai(j)
endif
endif
enddo
do j=1,n
fai(j)=dum(j)
enddo
enddo
do i=1,n
write(*,*)fai(i)
enddo
write(*,*)' '
end subroutine
subroutine inkaihou(fai,t,dt,x,dx,xmax)
implicit none
double precision fai0(10000),fai(10000),fai2(10000),dum(10000),t,dt,x,dx,xmax,pi,r,a(10000,10000)
integer n,m,i,j
pi=3.14159265359d0
x=0.0d0
n=dble(xmax/dx)
do i=1,n
x=x+dx
fai0(i)=sin(pi*x/2.0d0)
enddo
r=dt/(dx**2.0d0)
do i=1,n
a(i,i)=1+2*r
if(i<n)then
a(i,i+1)=-r
endif
if(i>1)then
a(i,i-1)=-r
endif
enddo
m=dble(t/dt)
do i=1,m
if(i==1) then
call gauss(a,fai,fai0,n)
else
call gauss(a,fai2,fai,n)
fai=fai2
endif
enddo
do i=1,n
write(*,*)fai(i)
enddo
write(*,*)' '
end subroutine
subroutine cranknicholson(fai,t,dt,x,dx,xmax)
implicit none
double precision fai0(10000),fai(10000),fai2(10000),dum(10000),t,dt,x
double precision dx,xmax,pi,r,a(10000,10000),a2(10000,10000),y(10000)
integer n,m,i,j
!pi=3.14159265359d0
!x=0.0d0
!n=dble(xmax/dx)
!do i=1,n
! x=x+dx
! fai0(i)=sin(pi*x/2)
!enddo
r=dt/(dx**2.0d0)
!a(1,1)=1+r
!write(*,*)r,n
!do i=1,n
! a(i,i)=1+r
! if(i<n)then
! a(i,i+1)=-r/2.0d0
! endif
! if(i>1)then
! a(i,i-1)=-r/2.0d0
! endif
!enddo
!do i=1,n
! a2(i,i)=1-r
! if(i<n)then
! a2(i,i+1)=r/2.0d0
! endif
! if(i>1)then
! a2(i,i-1)=r/2.0d0
! endif
!enddo
m=dble(t/dt)
!do i=1,m
!if(i==1) then
!call ax(a2,fai0,n,y) !{n*nの行列A}*{1*nの行列x}
!call gauss(a,fai,y,n)
!else
!call ax(a2,fai,n,y)
!call gauss(a,fai2,y,n)
!fai=fai2
!endif
!enddo
!do i=1,n
!write(*,*)fai(i)
!enddo
write(*,*)' '
end subroutine
subroutine ax(a,x,n,ans)!{n*nの行列A}*{1*nの行列x}
implicit none
double precision a(10000,10000),x(10000),dum(10000),ans(10000)
integer i,j,n
dum=0
do i=1,n
do j=1,n
dum(i)=dum(i)+x(j)*a(i,j)
enddo
enddo
do i=1,n
ans(i)=dum(i)
enddo
end subroutine
subroutine gauss(a,x,y,n)
implicit none
double precision a(10000,10000),x(10000),y(10000),kari,kari2,dum(10000,10000)
integer n,i,j,k,h,p
do i=1,n
do j=1,n
dum(i,j)=a(i,j)
enddo
enddo
h=0
!前進消去43~70
do i=1,n
do j=i,n
if(dum(j,i)/=0) then
exit
endif
if(j==n) then
write(*,*)"Unique solution does not exist"
h=1
endif
enddo
enddo
if(h==0)then
do i=1,n
call pivotselection(i,dum,y,n)
do p=1,n
enddo
do j=i,n
y(i)=y(i)/dum(i,i)
kari=dum(i,i)
if(j<n)then
kari2=dum(j+1,i)
endif
do k=i,n
dum(j,k)=dum(j,k)/kari
if(j<n)then
dum(j+1,k)=dum(j+1,k)-kari2*dum(i,k)
if(k==i)then
y(j+1)=y(j+1)-kari2*y(i)
endif
endif
enddo
enddo
enddo
endif
!後退代入
x(n)=y(n)/dum(n,n)
do i=1,n-1
x(n-i)=y(n-i)/dum(n-i,n-i)
do j=1,i
x(n-i)=x(n-i)-dum(n-i,n-j+1)*x(n-j+1)/dum(n-i,n-i)
enddo
enddo
end subroutine
subroutine pivotselection(s,a,y,n)
implicit none
double precision a(10000,10000),y(10000),saidai,kari
integer i,j,s,n!jは最大の係数をもつのは何行目かを保存、sは対角線上の成分の添え字、nはn*n行列のn
saidai=abs(a(s,s))
j=s
do i=s,n-1
if(saidai<abs(a(i+1,s))) then
saidai=a(i+1,s)
j=i+1
endif
enddo
do i=s,n
kari=a(s,i)
a(s,i)=a(j,i)
a(j,i)=kari
enddo
kari=y(s)
y(s)=y(j)
y(j)=kari
end subroutine
コードを貼るの忘れてました
[code]
program main
implicit none
double precision fai(10000),t,dt,x,dx,xmax
dx=0.1d0
xmax=2.0d0
dt=0.001d0
t=0.05d0
call youkaihou(fai,t,dt,x,dx,xmax)
call inkaihou(fai,t,dt,x,dx,xmax)
!call cranknicholson(fai,t,dt,x,dx,xmax)
end program main
subroutine youkaihou(fai,t,dt,x,dx,xmax)!t秒後の温度を求める。
implicit none
double precision fai0(10000),fai(10000),dum(10000),t,dt,x,dx,xmax,pi,r
integer n,m,i,j
pi=3.14159265359d0
x=0.0d0
n=dble(xmax/dx)
do i=1,n
x=x+dx
fai0(i)=sin(pi*x/2.0d0)
enddo
m=dble(t/dt)
r=dt/(dx**2.0d0)
do i=1,m
do j=1,n
if(i==1)then
if(j==1) then
dum(j)=(1.0d0-2.0d0*r)*fai0(j)+r*fai0(j+1)
else if(j==n) then
dum(j)=r*fai0(n-1)+(1.0d0-2.0d0*r)*fai0(n)
else!1<j<n-1のとき
dum(j)=r*fai0(j-1)+(1.0d0-2.0d0*r)*fai0(j)+r*fai0(j)
endif
else !i>2のとき
if(j==1) then
dum(j)=(1.0d0-2.0d0*r)*fai(j)+r*fai(j+1)
else if(j==n) then
dum(j)=r*fai(n-1)+(1.0d0-2.0d0*r)*fai(n)
else!1<j<n-1のとき
dum(j)=r*fai(j-1)+(1.0d0-2.0d0*r)*fai(j)+r*fai(j)
endif
endif
enddo
do j=1,n
fai(j)=dum(j)
enddo
enddo
do i=1,n
write(*,*)fai(i)
enddo
write(*,*)' '
end subroutine
subroutine inkaihou(fai,t,dt,x,dx,xmax)
implicit none
double precision fai0(10000),fai(10000),fai2(10000),dum(10000),t,dt,x,dx,xmax,pi,r,a(10000,10000)
integer n,m,i,j
pi=3.14159265359d0
x=0.0d0
n=dble(xmax/dx)
do i=1,n
x=x+dx
fai0(i)=sin(pi*x/2.0d0)
enddo
r=dt/(dx**2.0d0)
do i=1,n
a(i,i)=1+2*r
if(i<n)then
a(i,i+1)=-r
endif
if(i>1)then
a(i,i-1)=-r
endif
enddo
m=dble(t/dt)
do i=1,m
if(i==1) then
call gauss(a,fai,fai0,n)
else
call gauss(a,fai2,fai,n)
fai=fai2
endif
enddo
do i=1,n
write(*,*)fai(i)
enddo
write(*,*)' '
end subroutine
subroutine cranknicholson(fai,t,dt,x,dx,xmax)
implicit none
double precision fai0(10000),fai(10000),fai2(10000),dum(10000),t,dt,x
double precision dx,xmax,pi,r,a(10000,10000),a2(10000,10000),y(10000)
integer n,m,i,j
!pi=3.14159265359d0
!x=0.0d0
!n=dble(xmax/dx)
!do i=1,n
! x=x+dx
! fai0(i)=sin(pi*x/2)
!enddo
r=dt/(dx**2.0d0)
!a(1,1)=1+r
!write(*,*)r,n
!do i=1,n
! a(i,i)=1+r
! if(i<n)then
! a(i,i+1)=-r/2.0d0
! endif
! if(i>1)then
! a(i,i-1)=-r/2.0d0
! endif
!enddo
!do i=1,n
! a2(i,i)=1-r
! if(i<n)then
! a2(i,i+1)=r/2.0d0
! endif
! if(i>1)then
! a2(i,i-1)=r/2.0d0
! endif
!enddo
m=dble(t/dt)
!do i=1,m
!if(i==1) then
!call ax(a2,fai0,n,y) !{n*nの行列A}*{1*nの行列x}
!call gauss(a,fai,y,n)
!else
!call ax(a2,fai,n,y)
!call gauss(a,fai2,y,n)
!fai=fai2
!endif
!enddo
!do i=1,n
!write(*,*)fai(i)
!enddo
write(*,*)' '
end subroutine
subroutine ax(a,x,n,ans)!{n*nの行列A}*{1*nの行列x}
implicit none
double precision a(10000,10000),x(10000),dum(10000),ans(10000)
integer i,j,n
dum=0
do i=1,n
do j=1,n
dum(i)=dum(i)+x(j)*a(i,j)
enddo
enddo
do i=1,n
ans(i)=dum(i)
enddo
end subroutine
subroutine gauss(a,x,y,n)
implicit none
double precision a(10000,10000),x(10000),y(10000),kari,kari2,dum(10000,10000)
integer n,i,j,k,h,p
do i=1,n
do j=1,n
dum(i,j)=a(i,j)
enddo
enddo
h=0
!前進消去43~70
do i=1,n
do j=i,n
if(dum(j,i)/=0) then
exit
endif
if(j==n) then
write(*,*)"Unique solution does not exist"
h=1
endif
enddo
enddo
if(h==0)then
do i=1,n
call pivotselection(i,dum,y,n)
do p=1,n
enddo
do j=i,n
y(i)=y(i)/dum(i,i)
kari=dum(i,i)
if(j<n)then
kari2=dum(j+1,i)
endif
do k=i,n
dum(j,k)=dum(j,k)/kari
if(j<n)then
dum(j+1,k)=dum(j+1,k)-kari2*dum(i,k)
if(k==i)then
y(j+1)=y(j+1)-kari2*y(i)
endif
endif
enddo
enddo
enddo
endif
!後退代入
x(n)=y(n)/dum(n,n)
do i=1,n-1
x(n-i)=y(n-i)/dum(n-i,n-i)
do j=1,i
x(n-i)=x(n-i)-dum(n-i,n-j+1)*x(n-j+1)/dum(n-i,n-i)
enddo
enddo
end subroutine
subroutine pivotselection(s,a,y,n)
implicit none
double precision a(10000,10000),y(10000),saidai,kari
integer i,j,s,n!jは最大の係数をもつのは何行目かを保存、sは対角線上の成分の添え字、nはn*n行列のn
saidai=abs(a(s,s))
j=s
do i=s,n-1
if(saidai<abs(a(i+1,s))) then
saidai=a(i+1,s)
j=i+1
endif
enddo
do i=s,n
kari=a(s,i)
a(s,i)=a(j,i)
a(j,i)=kari
enddo
kari=y(s)
y(s)=y(j)
y(j)=kari
end subroutine
[/code]