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

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

zengji630

fcode 发表于 2014-12-13 09:25
只要这条语句执行了，赋值语句一定会成功。

谢谢你的解答

fcode

只要这条语句执行了，赋值语句一定会成功。

zengji630

vvt 发表于 2014-12-11 23:30
我说寻找 0.83083 的小数形式，并不是说，改成小数形式就可以解决问题的。
而是可以看出它是一个偶函数（只 ...

嗯嗯，知道啦。谢谢的解答
我想请问一下
SCO2  = 0.063*PaCO20**0.83083 /(1+0.63*PaCO20**0.83083+0.275*ST(PaO20+COMM*PaCO0) )
CaCO2 = AbCO2*PaCO20 + XNCO2*SCO2
CICO2(J) = CaCO2
CVCO2=CICO2(J)这里CVCO2有没有被赋值成功？

vvt

我说寻找 0.83083 的小数形式，并不是说，改成小数形式就可以解决问题的。
而是可以看出它是一个偶函数（只考虑实数域），所以可以在 x 上取一个 abs(x) 绝对值。

比如：

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

SCO2  = 0.063*abs(PaCO20)**(50000./60181)
     +  /(1+0.63*abs(PaCO20)**(50000./60181)+0.275*ST(PaO20+COMM*PaCO0) )

zengji630

zengji630 发表于 2014-12-11 10:05
谢谢你的帮助。我检查了一下dFunCO2,也重新调整了一下函数的计算顺序，还是出现一样的结果。是不是我的迭 ...

代码和附件都更新了一下

zengji630

vvt 发表于 2014-12-10 20:02
i , j , k 循环是。当 i = 1 , j = 2 , k = 2 时，
PaCO2 = PaCO20 - FunCO2(PaCO20)/ dFunCO2(PaCO20)
结 ...

谢谢你的帮助。我检查了一下dFunCO2,也重新调整了一下函数的计算顺序，还是出现一样的结果。是不是我的迭代有问题呢？

vvt

i , j , k 循环是。当 i = 1 , j = 2 , k = 2 时，
PaCO2 = PaCO20 - FunCO2(PaCO20)/ dFunCO2(PaCO20)
结果大概是 -2.0 导致负数

在上一次迭代时，PaCO20 = 55.41635，FunC02(PaCO20) = 342.9697 , dFunCO2 = 5.932721

所以 55.41635 - 342.9697 / 5.932721 = -2

zengji630

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

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

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

vvt

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

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

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

没有别的办法。

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

zengji630

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

FunCO2 函数中

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