大家好，负数取小数指数的问题

zengji630

[Fortran] 纯文本查看 复制代码

        Program Gas exchange
C        real AaO2, AbO2, XNO2, DLO2, O2MM, PIO2,PaCO,CVO2
C        real AaCO,AbCO, XNCO, DLCO, COMM, PICO, PaO2,CVCO
C        real AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2

        COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2

        COMMON /parametersCO/ AaCO,AbCO, XNCO, DLCO, COMM, PICO,PaO20,CVCO
        
        COMMON /parametersCO2/ AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2

      COMMON /VAR/  QV, QB, QE, QI ,Cos
        
      DIMENSION PBloodO2(20)
        DIMENSION PBloodCO(20)
        DIMENSION PBloodCO2(20)

        DIMENSION CIO2(20)
        DIMENSION CICO(20)
        DIMENSION CICO2(20)

!        Read the constants

        OPEN  (UNIT=1, FILE='CONSTANT.dat', STATUS='OLD')
        OPEN  (UNIT=2, FILE='PAlveolarO2.dat', STATUS='UNKNOWN') 
      OPEN  (UNIT=3, FILE='PAlveolarCO.dat', STATUS='UNKNOWN')
        OPEN  (UNIT=4, FILE='PAlveolarCO2.dat', STATUS='UNKNOWN')
        OPEN  (UNIT=5, FILE='COHB.dat', STATUS='UNKNOWN')
        OPEN  (UNIT=6, FILE='PBLOODO2.dat', STATUS='UNKNOWN') 
      OPEN  (UNIT=7, FILE='PBLOODCO.dat', STATUS='UNKNOWN')
        OPEN  (UNIT=8, FILE='PBLOODCO2.dat', STATUS='UNKNOWN')
        OPEN  (UNIT=9, FILE='CIO2.dat', STATUS='UNKNOWN')
      OPEN  (UNIT=10, FILE='CICO.dat', STATUS='UNKNOWN')
        OPEN  (UNIT=11, FILE='CICO2.dat', STATUS='UNKNOWN')
        
      write (2,*) "O2 partial pressure in Alveolar"
        write (3,*) "CO partial pressure in Alveolar"
        write (4,*) "CO2 partial pressure in Alveolar"
        write (5,*) "COHB percentage in blood"

         Read  (1,*)  AaO2,  AbO2,  XNO2,  DLO2,  O2MM
        Read  (1,*)  AaCO,  AbCO,  XNCO,  DLCO,  COMM  
        Read  (1,*)  AaCO2, AbCO2, XNCO2, DLCO2, CO2MM

!        Control parameters
        Write (*,*) "O2 exposure concentration =  (unit in mmHg)"
        Read  (*,*) PIO2
        Write (*,*) "CO exposure concentration =  (unit in mmHg)"
        Read  (*,*) PICO
        Write (*,*) "CO2 exposure concentration = (unit in mmHg)"
        Read  (*,*) PICO2
        Write (*,*) "Volume of blood in the body = (unit in ml)"
        Read  (*,*) QV
        Write (*,*) "Blood flow rate = (unit in ml/min)"
        Read  (*,*) QB
        Write (*,*) "Ventilation rate = (unit in ml/min)"
        Read  (*,*) QE

!        Control parameters
        
        Write (*,*) "Cycle number = "
        Read (*,*) ICycle
        
        Write (*,*) "Iterative number = "                            
        Read (*,*) Iteration

        Write (*,*) "Vascular Compartments = "
        Read (*,*) Cos

        CVO2 = 0.155742
        CVCO = 0.0000000001
        CVCO2 = 0.54382

        DO I = 1, ICycle
         
                DO J = 1, Cos
                                  
                        PaO20  = 0.0000000001
                         PaCO0  = 0.0000000001
                        PaCO20 = 0.0000000001

                        DO K = 1, Iteration
                           
                                PaO2  = PaO20  - FunO2(PaO20)  / dFunO2(PaO20)

                                PaCO  = PaCO0  - FunCO(PaCO0)  / dFunCO(PaCO0)

                                PaCO2 = PaCO20 - FunCO2(PaCO20)/ dFunCO2(PaCO20)

                        PaO20  = PaO2
                                PaCO0  = PaCO
                                PaCO20 = PaCO2

                        ENDDO

                        SO2  = ST(PaO20+COMM*PaCO0) / ( 1+(COMM*PaCO0)/PaO20 )
                        CaO2 = AbO2*PaO20 + XNO2*SO2
                        CIO2(J) = CaO2
                        CVO2=CIO2(J)

                        SCO  = ST(PaO20+COMM*PaCO0) / ( 1+PaO20/(COMM*PaCO0) ) 
                        CaCO = AbCO*PaCO0 + XNCO*SCO
                        CICO(J) = CaCO
                        CVCO=CICO(J)

                        SCO2  = 0.063*PaCO20**(5000/60181)
        +            /(1+0.63*PaCO20**(5000/60181)+0.275*ST(PaO20+COMM*PaCO0) )
                        CaCO2 = AbCO2*PaCO20 + XNCO2*SCO2
                          CICO2(J) = CaCO2
                        CVCO2=CICO2(J)

                        PBLOODO2(J)        = PaO20
                        PBLOODCO(J) = PaCO0 
                        PBLOODCO2(J) = PaCO20
        
              WRITE (6,*) J, PBLOODO2(J)
                        WRITE (7,*) J, PBLOODCO(J)                
                        WRITE (8,*) J, PBLOODCO2(J)

                WRITE (9,*) J, CIO2(J)
                        WRITE (10,*) J, CICO(J)        
                        WRITE (11,*) J, CICO2(J)

                END DO

                QI = QE + QB*( (CaO2-CVO2)+(CaCO-CVCO)
     +                  + (CaCO2-CVCO2) ) / 0.826184
                  CVO2 = CaO2 - 250/QB
                CVCO = CaCO + 0.007/QB
                CVCO2 = CaCO2 + 250/QB

                PaO2AVER=SUM(PBLOODO2)/Cos  
                PAlveolarO2 = ( QI*AaO2*PIO2 + DLO2*PaO2AVER) /(QE*AaO2 +DLO2)

            PaCOAVER=SUM(PBLOODCO)/Cos
                PAlveolarCO = ( QI*AaCO*PICO + DLCO*PaCOAVER )/(QE*AaCO +DLCO)

                PaCO2AVER=SUM(PBLOODCO2)/Cos 
            PAlveolarCO2 = ( QI*AaCO2*PICO2 + DLCO2*PaCO2AVER )
     +                                    / ( QE*AaCO2 + DLCO2 )

                Write (2,*) I, PAlveolarO2
                Write (3,*) I, PAlveolarCO 
                Write (4,*) I, PAlveolarCO2
                Write (5,*) I, SCO

        ENDDO

        End Program

        FUNCTION FunO2(x)

        COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2

        COMMON /VAR/  QV,QB, QE, QI,Cos

        FunO2 = ( QB*AbO2 + DLO2/Cos*QE*AaO2 / (DLO2/Cos+QE*AaO2) )*x
     +     +  XNO2*x*QB / (X+O2MM*PaCO0)
     +     *  ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) )
     +     /  ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )
     +     -  DLO2/Cos*QI*AaO2*PIO2 / (DLO2/Cos+QE*AaO2) - QB*CVO2

        END FUNCTION

        FUNCTION dFunO2(x)

        COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2

      COMMON /VAR/  QV,QB, QE, QI,Cos

        dFunO2 =  ( QB*AbO2 + DLO2/Cos*QE*AaO2 / (DLO2/Cos+QE*AaO2) )
     +     +  XNO2*QB
     +     *  ( ( (X+O2MM*PaCO0)**2 + 150 )
     +        * ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )
     +        + 2*X*(X+O2MM*PaCO0)
     +        * ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )
     +        -X*( (X+O2MM*PaCO0)**2 + 150 )
     +        *( 3*(X+O2MM*PaCO0)**2 + 150) )
     +     / ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )**2

        END FUNCTION
 
        FUNCTION FunCO(x)

        COMMON /parametersCO/ AaCO, AbCO, XNCO, DLCO, COMM,PICO,PaO20,CVCO
        
        COMMON /VAR/  QV, QB, QE, QI,Cos

        FunCO = ( QB*AbCO + DLCO/Cos*QE*AaCO / (DLCO/Cos+QE*AaCO) )*x
     +     +  XNCO*COMM*QB*x / (PaO20+COMM*x)
     +     *  ( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) )
     +     /  ( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) +23400 )
     +     -  DLCO/Cos*QI*AaCO*PICO  / (DLCO/Cos+QE*AaCO)        
     +     -  QB*CVCO

        END FUNCTION

        FUNCTION dFunCO(x)

        COMMON /parametersCO/ AaCO, AbCO, XNCO, DLCO, COMM,PICO,PaO20,CVCO
        
        COMMON /VAR/  QV,QB, QE, QI,Cos

        dFunCO = ( QB*AbCO + DLCO/Cos*QE*AaCO / (DLCO/Cos+QE*AaCO) )
     +     + XNCO*COMM*QB
     +     *( ( (PaO20+COMM*x)**2 + 2*COMM*x*(PaO20+COMM*x)+ 150)
     +     *  ( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) + 23400 )
     +         - ( (PaO20+COMM*x)**2 + 150 )
     +         *x*COMM*( 3*(PaO20+COMM*x)**2 + 150 ) )
     +     /  ( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) + 23400 )**2

        END FUNCTION

        FUNCTION FunCO2(X)

        COMMON /parametersCO2/ AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2

      COMMON /parametersCO/ AaCO,AbCO, XNCO, DLCO, COMM, PICO,PaO20,CVCO

        COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2

        COMMON /VAR/  QV,QB, QE, QI,Cos
        

        FunCO2 = ( QB*AbCO2 + (DLCO2/Cos*QE*AaCO2)/(DLCO2/Cos+QE*AaCO2))*X
        +      + QB*XNCO2*(0.063*x**(50000/60181))
     +          /( 1 + 0.63*x**(50000/60181) + 0.275*ST(PaO20+CO2MM*PaCO0) )
     +      - DLCO2/Cos*QI*AaCO2*PICO2 / (DLCO2/Cos+QE*AaCO2) 
     +      -        QB*CVCO2

        END FUNCTION

        FUNCTION dFunCO2(X)

        COMMON /parametersCO2/ AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2

      COMMON /parametersCO/ AaCO,AbCO, XNCO, DLCO, COMM, PICO,PaO20,CVCO

        COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2

        COMMON /VAR/  QV,QB, QE, QI,Cos
        

        dFunCO2 = QB*AbCO2 + (DLCO2/Cos*QE*AaCO2) / (DLCO2/Cos + QE*AaCO2)
          +           + QB*XNCO2*(0.05234*( 1 + 0.063*x**(50000/60181)
     +       + 0.275*ST(PaO20+CO2MM*PaCO0 ) - 0.063*x**(50000/60181) )
     +            / ( x**0.16917*( 1 + 0.063*x**(50000/60181)
     +        + 0.275*ST(PaO20+CO2MM*PaCO0) )**2 ) )

        END FUNCTION

        FUNCTION ST(x)

        ST = ( x**3 + 150*x ) / ( x**3 + 150*x + 23400 )        

        END FUNCTION

珊瑚虫

由于是固定格式，程序在我这边编译的时候，只有部分行长度超过了允许长度，需要往左侧挪一下，挪完之后，编译再没出现任何错误，至于运行期错误，楼主没有提供数据文件，所以无法知道是什么错误。
ps 本论坛的代码可以高亮显示的，楼主可以编辑下帖子 代码在 编辑菜单的 “<> ”输入

zengji630

哦，知道啦，谢谢你的指导。还有一个要读入的constant文件，屏幕输入的数据依次是0,0,149,5000,5000,5000,600,10。

我这里也能运行程序，就是只能循环两次自动跳出来，然后就出现First-chance exception in GAS.exe (KERNELBASE.DLL): 0xC0000090: Float Invalid Operation.检查了好几天了也不知道哪里出问题。希望你能帮我解答一下

vvt

9个屏幕输入，你只给了8个。

zengji630

vvt 发表于 2014-12-9 11:59
9个屏幕输入，你只给了8个。

0，0,149,5000,5000,5000,600,100,10谢谢

zengji630

珊瑚虫 发表于 2014-12-8 20:10
由于是固定格式，程序在我这边编译的时候，只有部分行长度超过了允许长度，需要往左侧挪一下，挪完之后，编 ...

哦，知道啦，谢谢你的指导。还有一个要读入的constant文件，屏幕输入的数据依次是0,0,149,5000,5000,5000,600,10。

我这里也能运行程序，就是只能循环两次自动跳出来，然后就出现First-chance exception in GAS.exe (KERNELBASE.DLL): 0xC0000090: Float Invalid Operation.检查了好几天了也不知道哪里出问题。希望你能帮我解答一下

vvt

你的问题在于：

FunCO2 函数中

[Fortran] 纯文本查看 复制代码

FunCO2 = ( QB*AbCO2 + (DLCO2/Cos*QE*AaCO2)/(DLCO2/Cos+QE*AaCO2))*x
     +      + QB*XNCO2*(0.063*x**0.83083)

这里有 x ** 0.83083 这个语句。 Fortran 中的指数为小数时，x 不能为负数。

这个 0.83083 是什么含义？可否写成分数的形式？

另外，为了安全，我建议你把 **2.0 都改为 **2，**3.0 都改为 **3
指数如果可以是整数，就不要写成小数。

zengji630

vvt 发表于 2014-12-9 23:20
你的问题在于：

FunCO2 函数中

我按照要求改了一下，还是出现一样的问题。0.83083是该函数中给定的值，没有整数形式，

vvt

负数不能开偶次方，所以负数也就不能取偶数倒数的指数。

这是数学问题决定的。不是 fortran 决定的。

解决方法：
1.寻求 0.83083 的分数形式。比如 多少除以多少这种。
2.保证 x 大于 0。

没有别的办法。

如果这两个办法你都无法做到，那我就帮不了你了。

zengji630

vvt 发表于 2014-12-10 16:04
负数不能开偶次方，所以负数也就不能取偶数倒数的指数。

这是数学问题决定的。不是 fortran 决定的。

非常感谢你的指导。第一点我改正了。第二点由于X是通过牛顿迭代法求出来的，不知道是不是负值。可是根据实际情况不可能是负的。