求高手看一下fortran程序？？？谢谢

xiaomaoxiaomao

小弟编写了一个PR状态方程的程序，用牛顿迭代法算v，运行没有问题，就是算的结果不如人意，求高手帮忙看一下，程序如下：

program PREOS
implicit none
  integer i,n
  real,parameter::R=83.14
  real,parameter::Tc=187.564
  real,parameter:c=45.592
  real,parameter::error=0.0000001
  real::p=319.144
  real::T=187.564
  real::w=0.01142
  real::m,b,a,Tr,f,f1,Pr,Pa
  real::v,v0
  real::z1,z2
  read(*,*)v0,n
  
  m=0.37464+1.54226*w-0.06992*w**2
  Tr=T/Tc
  a=0.45724*(R**2)*(Tc**2)*(1.0+m*(1.0-Tr**0.5))**2/Pc
  b=0.07780*R*Tc/Pc
  i=1
  do while(.true.)
   
     f=v0**3-(R*T/P-b)*v0**2+(a/P-2*b*R*T/P-3*b**2)*v0-b*(a/P-b*R*T/P-b**2)
     f1=3*v0**2-2*v0*(R*T/P-b)+(a/P-2*b*R*T/P-3*b**2)
     v=v0-f/f1
  if(abs(v-v0)<error)exit
  v0=v
  i=i+1
  if(i>=n)exit

   end do
   
   Pr=R*T/(V-b)
   Pa=-a/(v*(v+b)+b*(v-b))
   z1=Pr/Pc
   z2=Pa/Pc
   
   write(*,"('i=',i4)")i
   write(*,"('v=',f12.8)")v
   write(*,"('Pr=',f20.8)")Pr
   write(*,"('Pa=',f20.8)")Pa
   write(*,"('z1=',f20.8)")z1
   write(*,"('z2=',f20.8)")z2

   stop
end




其中PR状态方程如下所示：

xiaomaoxiaomao

1# xiaomaoxiaomao



第六行为
real,parameter::pc=45.592

guojunhang

你这代码属于不能拷贝的那一种，想帮你也用不上力啊

gdyu_yu

m=0.37464+1.54226*w-0.26992*w**2

xiaomaoxiaomao

4# gdyu_yu


那牛顿迭代应该没有问题吧？

xiaomaoxiaomao

3# guojunhang


可以拷贝吧，恩，我上边就是复制粘贴上去的，可以复制粘贴上来的呀？

铁道科学

看了一下，因为不大明白你的算法，暂时只能帮你把代码整理下，方便后面帮你的人看。

1. program PREOS
   
2. implicit none
   
3.   integer i,n
   
4.   real,parameter :: R=83.14
   
5.   real,parameter :: Tc=187.564
   
6.   real,parameter :: Pc=45.592
   
7.   real,parameter :: error=0.0000001
   
8.   real :: p=319.144
   
9.   real :: T=187.564
   
10.   real :: w=0.01142
    
11.   real :: m,b,a,Tr,f,f1,Pr,Pa
    
12.   real :: v,v0
    
13.   real :: z1,z2
    
14. 
    
15.   read(*,*) v0,n
    
16.   
    
17.   m  = 0.37464+1.54226*w-0.06992*w**2
    
18.   Tr = T/Tc
    
19.   a  = 0.45724*(R**2)*(Tc**2)*(1.0+m*(1.0-Tr**0.5))**2/Pc
    
20.   b  = 0.07780*R*Tc/Pc
    
21. 
    
22.   do i = 1,n
    
23.     f  = v0**3-(R*T/P-b)*v0**2+(a/P-2*b*R*T/P-3*b**2)*v0-b*(a/P-b*R*T/P-b**2)
    
24.     f1 = 3*v0**2-2*v0*(R*T/P-b)+(a/P-2*b*R*T/P-3*b**2)
    
25.     v  = v0-f/f1
    
26.     if(abs(v-v0)<error) exit
    
27.     v0 = v
    
28.   end do
    
29.    
    
30.   Pr = R*T/(V-b)
    
31.   Pa = -a/(v*(v+b)+b*(v-b))
    
32.   z1 = Pr/Pc
    
33.   z2 = Pa/Pc
    
34.    
    
35.   write(*,"('i=',i4)")i
    
36.   write(*,"('v=',f12.8)")v
    
37.   write(*,"('Pr=',f20.8)")Pr
    
38.   write(*,"('Pa=',f20.8)")Pa
    
39.   write(*,"('z1=',f20.8)")z1
    
40.   write(*,"('z2=',f20.8)")z2
    
41. 
    
42. end

复制代码

铁道科学

不知道下面的架构能不能满足你的要求，仅仅是一个大概的框架，未验证，请大家指教。

1. function Fx(x)
   
2.     函数体
   
3. end function

复制代码

主程序：

1. Program main
   
2.   implicit none
   
3.   real x,y,dx,yf,yt,nd
   
4. 
   
5.   nd = 1.0
   
6.   x  = 3.0
   
7.   yt = 9.0
   
8.   dx = 1.0
   
9.   yf = Fx(x)
   
10.   do while(dx>1e-6)
    
11.       x  = x+nd*dx
    
12.       y  = Fx(x)
    
13.       dx = nd*dx / (y-fy) * (yt-y)
    
14.   end do
    
15. 
    
16. end

复制代码

xiaomaoxiaomao

7# 铁道科学


谢谢。。。

xiaomaoxiaomao

8# 铁道科学

谢谢了。。

xiaomaoxiaomao

4# gdyu_yu

谢谢斑竹！！！！