24小时热门版块排行榜

返回列表

sshhyyy

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 59.1
帖子: 13
在线: 2.8小时
虫号: 2322361
注册: 2013-03-06
专业: 工程热力学

[求助] fortran计算多未知数非线性方程出错较少未知数时程序能通过

module inv_mat

contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine inv(A,invA,N)  !!计算逆矩阵
implicit real*8(a-z)
integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)
E=0
do i=1,n
   E(i,i)=1
end do
call mateq(A,E,invA,N,N)
end subroutine inv
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mateq(A,B,X,N,M) !!高斯列主元消去法计算矩阵方程
implicit real*8(a-z)
integer::N,M,i
real*8::A(N,N),B(N,M),X(N,M)
real*8::btemp(N),xtemp(N)
do i=1,M
   btemp=B(:,i)
   call elgauss(A,btemp,xtemp,N)
   X(:,i)=xtemp
end do

end subroutine mateq
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine elgauss(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)
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

   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 uptri(Aup,bup,x,n)
end subroutine elgauss
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine uptri(A,b,x,N) !!上三角方程组的回带方法
implicit real*8(a-z)
integer::N,j,i
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 uptri

end module inv_mat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module m_bfs

use inv_mat
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine solve(x0,N,tol,t_hsatoff,lbar,t_g4off,t_lsatoff,ggbar,lbarl,g_wh,g_wl,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)  !!方法函数
implicit real*8(a-z)
integer::i,n,itmax=50
real*8::x1(n),x0(n),y(n),f0(n),f1(n),dx(n)
real*8::v1(n),v2(n),v3(n)
real*8::H0(n,n),H1(n,n),df(n,n)
real*8::m1(n,n),m2(n,n),m3(n,n)

call jac(df,x0,t_hsatoff,dt_hjjoff,t_lsatoff,ggbar,lbar,lbarl,g_wh,g_wl,t_g4off,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)
call inv(df,H0,N)

!write(11,101)
!write(12,102)

do i=1,itmax
   call func(f0,x0,t_hsatoff,t_lsatoff,lbar,t_g4off,ggbar,lbarl,g_wh,g_wl,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)

   x1=x0-matmul(H0,f0)
   dx=x1-x0
   call func(f1,x1,t_hsatoff,t_lsatoff,lbar,t_g4off,ggbar,lbarl,g_wh,g_wl,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)
   y=f1-f0
   v1=matmul(H0,y)
   t1=vdot(y,v1,n)
   t2=vdot(dx,y,N)
   u=1d0+t1/t2

   m1=vvmat(dx,dx,n)*u
   m2=vvmat(dx,y,n)
   m2=matmul(m2,H0)
   v3=matmul(H0,y)
   m3=vvmat(v3,dx,n)
   t3=vdot(dx,y,n)
   H1=(m1-m2-m3)/t3
   H1=H0+H1
   write (12,103)i,H0
   x0=x1
   H0=H1
   write(11,104)i,x0
   write(*,*) x1(1)
   write(*,*) x0(1)
   dx2=dsqrt(dx(1)**2+dx(2)**2)
   if (dx2<tol) exit

end do
   write(*,*)x0(1)
101 format(/,t5,'bfs方法计算序列',/)
102 format(/,t5,'bfs方法h矩阵序列为:',/)
103 format(t5,'iter=',i4,/,<n>(<n>f16.10/),/)
104 format(i8,3f16.10)
end subroutine solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function vdot(a,b,N) !! 计算向量内积
integer::i,N
real*8::a(n),b(n),vdot
vdot=0
do i=1,N
   vdot=vdot+a(i)*b(i)
end do
end function vdot
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function vvmat(a,b,n)
integer::i,j,n
real*8::a(n),b(n),vvmat(n,n)
do i=1,n
   do j=1,n
   vvmat(i,j)=a(i)*b(j)
   end do
end do
end function vvmat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine func(f,x,t_hsatoff,t_lsatoff,lbar,t_g4off,ggbar,lbarl,g_wh,g_wl,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)  !!方程函数
implicit real*8(a-z)
real*8::x(14),f(14)

f(1)=-b1*ggbar*(t_g4off-x(1))+x(9)*(x(2)-t_hsatoff)
f(2)=-b2*exp(0.2*log(ggbar))*(t_g4off-x(1))+t_g4off-x(2)+x(1)-t_hsatoff
f(3)=-b3*ggbar*(t_g4off-x(1))+x(9)*(x(3)-x(4))
f(4)=-b4*exp(0.2*log(ggbar))*(t_g4off-x(1))+t_g4off-x(3)+x(1)-x(4)
f(5)=-x(9)*b5*(lbar+x(6))+ggbar*(x(1)-t_hsatoff-x(5))
f(6)=-exp(1/exp(0.2*log(ggbar))*log(b6))*x(5)+x(1)-t_hsatoff
f(7)=-b7*ggbar*(t_hsatoff+x(5)-x(7))+x(10)*(x(8)-t_lsatoff)
f(8)=-b8*exp(0.2*log(ggbar))*(t_hsatoff+x(5)-x(7))+t_hsatoff+x(5)-x(8)+x(7)-t_lsatoff
f(9)=b9*(t_hsatoff-x(6)-t_lsatoff)*x(9)+b9*x(9)*x(14)+ggbar*x(11)-ggbar*x(7)
f(10)=b10*(t_hsatoff-x(6)-t_lsatoff+x(14))*x(9)-exp(0.8*log(ggbar))*(x(7)-t_hsatoff+x(6)-t_lsatoff+x(14)+x(11))
f(11)=x(10)*b11*(lbarl+x(14))-ggbar*(x(11)-x(13)-t_lsatoff)
f(12)=exp(1/exp(0.2*log(ggbar))*log(b12))*x(13)-x(11)+t_lsatoff
f(13)=ggbar*(g_wh+g_wl)*(t_lsatoff+x(13)-x(12))-b13*(g_wh*x(9)+g_wl*x(10))*(t_lsatoff-x(14)-t_l1-273.15)
f(14)=exp(0.8*log(ggbar))*(g_wh+g_wl)*(x(13)+x(14)+x(12)-t_l1-273.15)-b14*(g_wh*x(9)+g_wl*x(10))*(t_lsatoff-x(14)-t_l1-273.15)

end subroutine func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine jac(df,x,t_hsatoff,dt_hjjoff,t_lsatoff,ggbar,lbar,lbarl,g_wh,g_wl,t_g4off,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)

implicit real*8(a-z)
real*8::x(14),df(14,14)

df(1,1)=b1*ggbar
df(1,2)=x(9)
df(1,3)=0d0
df(1,4)=0d0
df(1,5)=0d0
df(1,6)=0d0
df(1,7)=0d0
df(1,8)=0d0
df(1,9)=x(2)-t_hsatoff
df(1,10)=0d0
df(1,11)=0d0
df(1,12)=0d0
df(1,13)=0d0
df(1,14)=0d0

df(2,1)=b2*exp(0.2*log(ggbar))+1
df(2,2)=-1d0
df(2,3)=0d0
df(2,4)=0d0
df(2,5)=0d0
df(2,6)=0d0
df(2,7)=0d0
df(2,8)=0d0
df(2,9)=0d0
df(2,10)=0d0
df(2,11)=0d0
df(2,12)=0d0
df(2,13)=0d0
df(2,14)=0d0

df(3,1)=b3*ggbar
df(3,2)=0d0
df(3,3)=x(9)
df(3,4)=-x(9)
df(3,5)=0d0
df(3,6)=0d0
df(3,7)=0d0
df(3,8)=0d0
df(3,9)=x(3)-x(4)
df(3,10)=0d0
df(3,11)=0d0
df(3,12)=0d0
df(3,13)=0d0
df(3,14)=0d0

df(4,1)=b4*exp(0.2*log(ggbar))+1
df(4,2)=0d0
df(4,3)=-1d0
df(4,4)=-1d0
df(4,5)=0d0
df(4,6)=0d0
df(4,7)=0d0
df(4,8)=0d0
df(4,9)=0d0
df(4,10)=0d0
df(4,11)=0d0
df(4,12)=0d0
df(4,13)=0d0
df(4,14)=0d0

df(5,1)=ggbar
df(5,2)=0d0
df(5,3)=0d0
df(5,4)=0d0
df(5,5)=-ggbar
df(5,6)=-x(9)*b5
df(5,7)=0d0
df(5,8)=0d0
df(5,9)=-b5*(lbar+x(6))
df(5,10)=0d0
df(5,11)=0d0
df(5,12)=0d0
df(5,13)=0d0
df(5,14)=0d0

df(6,1)=1d0
df(6,2)=0d0
df(6,3)=0d0
df(6,4)=0d0
df(6,5)=-exp(1/exp(0.2*log(ggbar))*log(b6))
df(6,6)=0d0
df(6,7)=0d0
df(6,8)=0d0
df(6,9)=0d0
df(6,10)=0d0
df(6,11)=0d0
df(6,12)=0d0
df(6,13)=0d0
df(6,14)=0d0

df(7,1)=0d0
df(7,2)=0d0
df(7,3)=0d0
df(7,4)=0d0
df(7,5)=-b7*ggbar
df(7,6)=0d0
df(7,7)=b7*ggbar
df(7,8)=x(10)
df(7,9)=0d0
df(7,10)=x(8)-t_lsatoff
df(7,11)=0d0
df(7,12)=0d0
df(7,13)=0d0
df(7,14)=0d0

df(8,1)=0d0
df(8,2)=0d0
df(8,3)=0d0
df(8,4)=0d0
df(8,5)=-b8*exp(0.2*log(ggbar))+1
df(8,6)=0d0
df(8,7)=b8*exp(0.2*log(ggbar))+1
df(8,8)=-1d0
df(8,9)=0d0
df(8,10)=0d0
df(8,11)=0d0
df(8,12)=0d0
df(8,13)=0d0
df(8,14)=0d0

df(9,1)=0d0
df(9,2)=0d0
df(9,3)=0d0
df(9,4)=0d0
df(9,5)=0d0
df(9,6)=-b9*x(9)
df(9,7)=-ggbar
df(9,8)=0d0
df(9,9)=b9*(t_hsatoff-x(6)-t_lsatoff+x(14))
df(9,10)=0d0
df(9,11)=ggbar
df(9,12)=0d0
df(9,13)=0d0
df(9,14)=b9*x(9)

df(10,1)=0d0
df(10,2)=0d0
df(10,3)=0d0
df(10,4)=0d0
df(10,5)=0d0
df(10,6)=-b10*x(9)-exp(0.8*log(ggbar))
df(10,7)=-exp(0.8*log(ggbar))
df(10,8)=0d0
df(10,9)=b10*(t_hsatoff-x(6)-t_lsatoff)+b10*x(14)
df(10,10)=0d0
df(10,11)=-exp(0.8*log(ggbar))
df(10,12)=0d0
df(10,13)=0d0
df(10,14)=b10*x(9)-exp(0.8*log(ggbar))

df(11,1)=0d0
df(11,2)=0d0
df(11,3)=0d0
df(11,4)=0d0
df(11,5)=0d0
df(11,6)=0d0
df(11,7)=0d0
df(11,8)=0d0
df(11,9)=0d0
df(11,10)=b11*(lbarl+x(14))
df(11,11)=-ggbar
df(11,12)=0d0
df(11,13)=ggbar
df(11,14)=x(10)*b11

df(12,1)=0d0
df(12,2)=0d0
df(12,3)=0d0
df(12,4)=0d0
df(12,5)=0d0
df(12,6)=0d0
df(12,7)=0d0
df(12,8)=0d0
df(12,9)=0d0
df(12,10)=0d0
df(12,11)=-1d0
df(12,12)=0d0
df(12,13)=exp(1/exp(0.2*log(ggbar))*log(b12))
df(12,14)=0d0

df(13,1)=0d0
df(13,2)=0d0
df(13,3)=0d0
df(13,4)=0d0
df(13,5)=0d0
df(13,6)=0d0
df(13,7)=0d0
df(13,8)=0d0
df(13,9)=-b13*g_wh*(t_lsatoff-x(14)-t_l1-273.15)
df(13,10)=-b13*g_wl*(t_lsatoff-x(14)-t_l1-273.15)
df(13,11)=0d0
df(13,12)=-ggbar*(g_wh+g_wl)
df(13,13)=ggbar*(g_wh+g_wl)
df(13,14)=b13*(g_wh*x(9)+g_wl*x(10))

df(14,1)=0d0
df(14,2)=0d0
df(14,3)=0d0
df(14,4)=0d0
df(14,5)=0d0
df(14,6)=0d0
df(14,7)=0d0
df(14,8)=0d0
df(14,9)=-b14*g_wh*(t_lsatoff-x(14)-t_l1-273.15)
df(14,10)=-b14*g_wl*(t_lsatoff-x(14)-t_l1-273.15)
df(14,11)=0d0
df(14,12)=exp(0.8*log(ggbar))*(g_wh+g_wl)
df(14,13)=exp(0.8*log(ggbar))*(g_wh+g_wl)
df(14,14)=exp(0.8*log(ggbar))*(g_wh+g_wl)+b14*(g_wh*x(9)+g_wl*x(10))

end subroutine jac

end module m_bfs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program main

use inv_mat
use m_bfs


implicit real*8(a-z)
integer::N=14
real*8::x0(14)

!open(unit=11,file='result.txt')
!open(unit=12,file='hmatrix.txt')
x0=(/0.1d0,0.1d0,-0.1d0,0.1d0,0.1d0,0.1d0,0.1d0,0.1d0,-0.1d0,0.1d0,0.1d0,0.1d0,0.1d0,0.1d0/)
t_hsatoff=1
lbar=6
t_lsatoff=3
ggbar=4
t_g4off=5
lbarl=6
g_wh=7
g_wl=8
b1=2
b2=2
b3=2
b4=2
b5=2
b6=2
b7=2
b8=2
b9=9
b10=1
b11=3
b12=2
b13=4
b14=5
t_l1=6
call solve(x0,N,1d-8,t_hsatoff,lbar,t_g4off,t_lsatoff,ggbar,lbarl,g_wh,g_wl,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,t_l1)
write(*,*) x0(1)
write(*,*) x0(2)
write(*,*) x0(3)
read(*,*)
end program main

回复此楼