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