[Fortran] 纯文本查看 复制代码
program main
implicit none
real::a(3,3)=(/10,-1,-1,-1,10,-1,-2,-2,5/)
real::b(3)=(/7.2,8.3,4.2/)
real x(3)
real i,j,time,k,ii,o(3),keep(3),ans,ans1,ans2
call Jacobi(i,j,time,ans,o,keep,x,a,b)
call Gauss_Seidel(i,j,time,k,ii,o,keep,ans1,ans2,a,b,x)
pause
end program main
subroutine Jacobi(i,j,time,ans,o,keep,x,a,b)
implicit none
real ans,o(3),keep(3),x(3),a(3,3),b(3),i,j,time
x(1)=0
x(2)=0
x(3)=0
do time=1,10
do i=1,3
keep(i)=x(i)
ans=0
do j=1,3
if(j==i) cycle
ans=ans+a(i,j)*x(j)
end do
x(i)=(b(i)-ans)/a(i,i)
o(i)=abs(x(i)-keep(i))
end do
write(*,*) x(1),x(2),x(3)
if(o(1)<0.001 .and.o(2)<0.001 .and.o(3)<0.001 ) exit
write(*,*) 'This is the answer!'
end do
return
pause
end
subroutine Gauss_Seidel(i,j,time,k,ii,o,keep,ans1,ans2,a,b,x)
implicit none
real o(3),keep(3),a(3,3),b(3),x(3)
real ans1,ans2,i,j,time,k,ii
x(1:3)=0
do time=1,10
do i=1,3
keep(i)=x(i)
end do
ans1=0
ans2=0
do j=2,3
ans1=ans1+a(1,j)*x(j)
end do
x(1)=(b(1)-ans1)/a(1,1)
x(2)=(b(2)-a(2,1)*x(1)-a(2,3)*x(3))/a(2,2)
do k=1,2
ans2=ans2+a(3,k)*x(k)
end do
x(3)=(b(3)-ans2)/a(3,3)
do ii=1,3
o(ii)=abs(x(ii)-keep(ii))
end do
write(*,*) x(1),x(2),x(3)
if(o(1)<0.001 .and.o(2)<0.001 .and.o(3)<0.001) exit
write(*,*) 'This is the answer!'
end do
return
pause
end