Fortran Coder

查看: 3555|回复: 21

[数值问题] fortran 程序运行不守恒是什么问题?急

[复制链接]

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
发表于 2023-8-22 10:20:26 | 显示全部楼层 |阅读模式
大家好,在这里发一个求助帖。
问题是这样的:在数学上可以验证出来三个函数表达式相加(三个函数表达式相加的结果与任何数都无关,之和就是守恒的)是绝对守恒的,但是用Fortran程序运行出来的结果不守恒,请问一下大家可能是什么问题?

258

帖子

0

主题

0

精华

版主

World Analyser

F 币
713 元
贡献
508 点

新人勋章美女勋章元老勋章热心勋章规矩勋章管理勋章

QQ
发表于 2023-8-22 11:53:57 | 显示全部楼层
函数数值不稳定,浮点数误差等等。

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
 楼主| 发表于 2023-8-23 10:44:07 | 显示全部楼层
那请问这个问题要怎么解决啊,怎么才能提高精度啊?我在写程序的时候都已经用的是双精度数了,求解,

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
 楼主| 发表于 2023-8-23 10:45:41 | 显示全部楼层
kyra 发表于 2023-8-22 11:53
函数数值不稳定,浮点数误差等等。

或者说,根据经验,哪个编译器的运算精度更高?

258

帖子

0

主题

0

精华

版主

World Analyser

F 币
713 元
贡献
508 点

新人勋章美女勋章元老勋章热心勋章规矩勋章管理勋章

QQ
发表于 2023-8-23 14:50:52 | 显示全部楼层
编译器的差别几乎可以忽略,主要还是要针对算法本身做优化,这方面的经验不容易总结和归纳。
如果你的算法不是很复杂,可以给出来,大家有针对性的探讨。

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
 楼主| 发表于 2023-8-24 09:27:14 | 显示全部楼层
kyra 发表于 2023-8-23 14:50
编译器的差别几乎可以忽略,主要还是要针对算法本身做优化,这方面的经验不容易总结和归纳。
如果你的算法 ...

其实程序是很简单的,只有几个解析式,就是最基本的加减乘除和乘方运算,连微积分都不涉及。但是它运行结果就是不守恒了。
而且我用Excel表格计算,用其他的电脑试,包括用其他的程序语言进行实验,出现的结果相同,都是不守恒的。

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
 楼主| 发表于 2023-8-24 09:42:01 | 显示全部楼层
八爪鱼 发表于 2023-8-23 10:45
或者说,根据经验,哪个编译器的运算精度更高?

[Fortran] 纯文本查看 复制代码
program sasm          !model with P32
    implicit real*8(a-h,o-z)
    integer,parameter::n=15
    integer,parameter::i=10
    real*8,parameter::a=1d0
    real*8,dimension(n)::p12,p21,p23,p32
    real*8,dimension(n)::b,c
    real*8,dimension(n)::r1,r2
    real*8,dimension(n)::c1,c2,c3
    real*8,dimension(n)::m1,m2,m3
    real*8,dimension(n)::x1,x2,x3
    real*8,dimension(n)::yn1,yn2,yn3,yn
    
    open(1,file='D:/program/0809/frequency12.txt')
    open(2,file='D:/program/0809/frequency21.txt')
    open(3,file='D:/program/0809/frequency23.txt')
    open(5,file='D:/program/0809/frequency32.txt')
    open(7,file='D:/program/0809/1200.txt')
    
    do ix=1,n
        read(1,*)xtemp,ytemp,p12(ix)
        read(2,*)xtemp,ytemp,p21(ix)
        read(3,*)xtemp,ytemp,p23(ix)
        read(5,*)xtemp,ytemp,p32(ix)
    end do
    close(1)
    close(2)
    close(3)
    close(5)
    
    b(i)=p12(i)+p21(i)+p23(i)+p32(i)
    c(i)=p12(i)*p23(i)+p12(i)*p32(i)+p21(i)*p32(i)
    
    r1(i)=(-1d0*b(i)+dsqrt(b(i)**(2d0)-4d0*a*c(i)))/(2d0*a)
    r2(i)=(-1d0*b(i)-dsqrt(b(i)**(2d0)-4d0*a*c(i)))/(2d0*a)

    c1(i)=p12(i)*(r2(i)+p12(i)+p21(i))/(r1(i)*(r1(i)-r2(i)))
    c2(i)=p12(i)*(r1(i)+p12(i)+p21(i))/(r2(i)*(r2(i)-r1(i)))
    c3(i)=1d0-c1(i)-c2(i)
    
    m1(i)=(r1(i)+p12(i))/p21(i)
    m2(i)=(r2(i)+p12(i))/p21(i)
    m3(i)=p12(i)/p21(i)
    
    x1(i)=(r1(i)**(2d0)+r1(i)*(p12(i)+p21(i)+p23(i))+p12(i)*p23(i))/(p21(i)*p32(i))
    x2(i)=(r2(i)**(2d0)+r2(i)*(p12(i)+p21(i)+p23(i))+p12(i)*p23(i))/(p21(i)*p32(i))
    x3(i)=(p12(i)*p23(i))/(p21(i)*p32(i))
    
    
    do t=0d0,10000000d0,10000d0
        yn1(i)=c1(i)*dexp(r1(i)*t)+c2(i)*dexp(r2(i)*t)+c3(i)
        yn2(i)=c1(i)*m1(i)*dexp(r1(i)*t)+c2(i)*m2(i)*dexp(r2(i)*t)+c3(i)*m3(i)
        yn3(i)=c1(i)*x1(i)*dexp(r1(i)*t)+c2(i)*x2(i)*dexp(r2(i)*t)+c3(i)*x3(i)
        
        yn(i)=yn1(i)+yn2(i)+yn3(i)
        
        write(*,*)t,yn(i)
        write(7,*)yn1(i),yn2(i),yn3(i),yn(i)    
      end do
    stop
    end program sasm

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
 楼主| 发表于 2023-8-24 09:44:41 | 显示全部楼层
kyra 发表于 2023-8-23 14:50
编译器的差别几乎可以忽略,主要还是要针对算法本身做优化,这方面的经验不容易总结和归纳。
如果你的算法 ...

您好,我已经把我的代码贴出来了,里面可能有误差的有一个e指数的运算,能方便帮我看一下吗?程序整体上来说是比较简单的

213

帖子

2

主题

0

精华

宗师

F 币
2126 元
贡献
875 点

规矩勋章

发表于 2023-8-24 22:20:37 | 显示全部楼层
本帖最后由 风平老涡 于 2023-8-24 22:31 编辑
八爪鱼 发表于 2023-8-24 09:44
您好,我已经把我的代码贴出来了,里面可能有误差的有一个e指数的运算,能方便帮我看一下吗?程序整体上 ...

有几个问题:
1)循环变量建议用整数型
2)从31行到47行没有在循环内,所有的计算都没有意义。
3)从51行到58行数组引用与循环变量名不符,

14

帖子

1

主题

0

精华

入门

F 币
63 元
贡献
34 点
 楼主| 发表于 2023-8-25 09:25:33 | 显示全部楼层
风平老涡 发表于 2023-8-24 22:20
有几个问题:
1)循环变量建议用整数型
2)从31行到47行没有在循环内,所有的计算都没有意义。

您好,非常感谢您提的建议。
我在这里给您回复一下:
1、关于您提到的循环变量用整数型,循环变量t是要参与计算的,其他的都是浮点型,如果只有t是整型量的话,不知道对计算结果是否有影响,这个我再试一下;
2、31行到47行没有在循环内,这个是这样的,主要是31行到47行中有变量i,而变量i我在第4行已经定义了。这个程序需要提前读取文件,文件里面有15个不同的数。所以定义了一个变量i,让其从1到15,每次在给定i的情况下,再让程序进行下面的循环计算,i的值我是每次手动改的;
3、数组引用与循环变量名不符?可以再解释一下吗?我这边没懂您说的是什么意思。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-6-19 19:02

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表