Fortran Coder

标题: 大家好,负数取小数指数的问题 [打印本页]

作者: zengji630    时间: 2014-12-8 17:11
标题: 大家好,负数取小数指数的问题
[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


CONSTANT.dat

122 Bytes, 下载次数: 1


作者: 珊瑚虫    时间: 2014-12-8 20:10
由于是固定格式,程序在我这边编译的时候,只有部分行长度超过了允许长度,需要往左侧挪一下,挪完之后,编译再没出现任何错误,至于运行期错误,楼主没有提供数据文件,所以无法知道是什么错误。
ps 本论坛的代码可以高亮显示的,楼主可以编辑下帖子 代码在 编辑菜单的 “<> ”输入
作者: zengji630    时间: 2014-12-9 09:27
哦,知道啦,谢谢你的指导。还有一个要读入的constant文件,屏幕输入的数据依次是0,0,149,5000,5000,5000,600,10。

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

CONSTANT.dat

151 Bytes, 下载次数: 4


作者: vvt    时间: 2014-12-9 11:59
9个屏幕输入,你只给了8个。
作者: zengji630    时间: 2014-12-9 16:09
vvt 发表于 2014-12-9 11:59
9个屏幕输入,你只给了8个。

0,0,149,5000,5000,5000,600,100,10谢谢
作者: zengji630    时间: 2014-12-9 16:10
珊瑚虫 发表于 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    时间: 2014-12-9 23:20
你的问题在于:

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    时间: 2014-12-10 15:18
vvt 发表于 2014-12-9 23:20
你的问题在于:

FunCO2 函数中

我按照要求改了一下,还是出现一样的问题。0.83083是该函数中给定的值,没有整数形式,
作者: vvt    时间: 2014-12-10 16:04
本帖最后由 vvt 于 2014-12-10 16:06 编辑

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

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

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

没有别的办法。

如果这两个办法你都无法做到,那我就帮不了你了。
作者: zengji630    时间: 2014-12-10 17:30
vvt 发表于 2014-12-10 16:04
负数不能开偶次方,所以负数也就不能取偶数倒数的指数。

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

非常感谢你的指导。第一点我改正了。第二点由于X是通过牛顿迭代法求出来的,不知道是不是负值。可是根据实际情况不可能是负的。
作者: vvt    时间: 2014-12-10 20:02
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    时间: 2014-12-11 10:05
vvt 发表于 2014-12-10 20:02
i , j , k 循环是。当 i = 1 , j = 2 , k = 2 时,
PaCO2 = PaCO20 - FunCO2(PaCO20)/ dFunCO2(PaCO20)
结 ...

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

代码和附件都更新了一下
作者: vvt    时间: 2014-12-11 23:30
本帖最后由 vvt 于 2014-12-11 23:50 编辑

我说寻找 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    时间: 2014-12-12 16:05
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有没有被赋值成功?
作者: fcode    时间: 2014-12-13 09:25
只要这条语句执行了,赋值语句一定会成功。
作者: zengji630    时间: 2014-12-16 15:08
fcode 发表于 2014-12-13 09:25
只要这条语句执行了,赋值语句一定会成功。

谢谢你的解答




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2