module M_GAUSS
!高斯列主元消去法模块
contains
subroutine LINEQ(A,B,X,N)
!高斯列主元消去法
implicit real*8(A-Z)
integer::I,K,N
integer::ID_MAX !主元素标号
real*8::A(N,N),B(N),X(N)
real*8::AUP(N,N),BUP(N)
!A,B为增广矩阵
real*8::AB(N,N+1)
real*8::VTEMP1(N+1),VTEMP2(N+1)
AB(1:N,1:N)=A
AB(:,N+1)=B
!列主元消去法的核心部分
do K=1,N-1
ELMAX=DABS(AB(K,K))
ID_MAX=K
!这段程序不是为了赋值最大元素给ELMAX,而是为了找出最大元素对应的标号
do I=K+1,N
if(dabs(AB(I,K))>ELMAX) then
ELMAX=AB(I,K)
ID_MAX=I
end if
end do
!交换两行元素,其他不变
VTEMP1=AB(K,:)
VTEMP2=AB(ID_MAX,:)
AB(K,:)=VTEMP2
AB(ID_MAX,:)=VTEMP1
!以上交换两行元素,然后开始消元
do I=K+1,N
TEMP=AB(I,K)/AB(K,K)
AB(I,:)=AB(I,:)-TEMP*AB(K,:)
end do
end do
AUP(:,:)=AB(1:N,1:N)
BUP(:)=AB(:,N+1)
!调用上三角方程组的回代方法
call UP_TRI(AUP,BUP,X,N)
end subroutine LINEQ
subroutine UP_TRI(A,B,X,N)
!上三角方程组的回代方法
implicit real*8(A-Z)
integer::I,J,N
real*8::A(N,N),B(N),X(N)
X(N)=B(N)/A(N,N)
!回代部分
do I=N-1,1,-1
X(I)=B(I)
do J=I+1,N
X(I)=X(I)-A(I,J)*X(J)
end do
X(I)=X(I)/A(I,I)
end do
end subroutine UP_TRI
end module M_GAUSS
module M_NEWTON
!关于方程的模块
contains
subroutine SOLVE()
!计算非线性方程组
use M_GAUSS
implicit real*8(A-Z)
integer::I,ITMAX=50
integer::N=2
real*8::X(2),F(2),DX(2)
real*8::DF(2,2)
!ITMAX为最大允许迭代次数
!N为方程组维数
!DF为偏导数矩阵
open(unit=11,file="jieguo.txt")
write(11,100)
100 format(/,T6,"Newton法计算非线性方程组迭代序列:",/)
X=(/1.6,1.2/)
TOL=1D-9
do I=1,ITMAX
call FUNC(F,X)
call JAC(DF,X)
call LINEQ(DF,-F,DX,N)
X=X+DX
write(11,200)I,X
200 format(I5,2F16.10)
!判断计算精度,满足误差要求时推出循环
DX2=DSQRT(DX(1)**2+DX(2)**2)
if(DX2<TOL) exit
end do
end subroutine SOLVE
subroutine FUNC(F,X)
!方程函数
implicit real*8(A-Z)
real*8::X(2),F(2)
F(1) = X(1)**2 + X(2)**2 - 4
F(2) = X(1)**2 - X(2)**2 - 1
end subroutine FUNC
subroutine JAC(DF,X)
!偏导数矩阵
implicit real*8(A-Z)
real*8::X(2),DF(2,2)
DF(1,1) = 2*X(1)
DF(1,2) = 2*X(2)
DF(2,1) = 2*X(1)
DF(2,2 )= -2*X(2)
end subroutine JAC
end module M_NEWTON
program newton
use M_NEWTON
!调用Newton法中的函数
call SOLVE()
end program newton
评论5