Fortran Coder

标题: 这究竟是算法不对还是语法错了？求指点！ [打印本页]

作者: 武清河 时间: 2014-4-6 22:24
标题: 这究竟是算法不对还是语法错了？求指点！
本帖最后由武清河于 2014-4-7 10:36 编辑

我把算法用在数组中求一组解，程序也能够走，只是结果不对，请大家进来看看，指点指点，不胜感激！代码如下：

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





program main

real,dimension(8)::z1,v1,q1,p1

real::t=10800,z2_jiashe,q_pj,f,h,g,z2_suan

real,dimension(7)::z2,v2,q2,p2

open(1,file='1.dat')

open(2,file='2.out')

read(1,*)q1,q2

z1(1)=38.0

p1(1)=174.0

do i=1,8

   v1(i)=g(z1(i))

   q_pj=(q1(i)+q2(i))/2

       do j=1,10000

           z2_jiashe=z1(i)+0.001*j

           p2(i)=f(z2_jiashe)

           v2(i)=q_pj*t-(p1(i)+p2(i))/2*t+v1(i)

           z2_suan=h(v2(i))

                if(abs(z2_jiashe-z2_suan)<1e-3)then

                    z2(i)=z2_jiashe

                    v2(i)=g(z2(i))

                    z1(i+1)=z2(i)

                    p1(i+1)=p2(i)

                endif

       enddo

enddo

write(2, '(f30.9)'),z2

end program main

function f(x) result(f_re)

real::x,f_re

f_re=32.93*x**2-2344*x+41710

end function f

function g(x) result(g_re)

real::x,g_re

g_re=1.911e-4*x**7.297

end function g

function h(x) result(h_re)

real::x,h_re

h_re=3.24*x**0.1369

end function h

1.dat

83 Bytes, 下载次数: 3

输入文件

作者: fcode 时间: 2014-4-7 08:16
1.请给出输入文件（或其内容）。1.dat
2.建议：写代码要搞好排版，缩进。
3.目测你的代码错误很多，不知道你是如何执行出结果来的。
比如，
real,dimension(8)::z1,v1,q1,p1
real,dimension(7)::z2,v2,q2,p2
这些数组都是8，7的长度。
但是循环是
do i=1,11
后面的各种越界。

作者: 武清河 时间: 2014-4-7 10:37

fcode 发表于 2014-4-7 08:16
1.请给出输入文件（或其内容）。1.dat
2.建议：写代码要搞好排版，缩进。
3.目测你的代码错误很多，不知道 ...

已修改，初次发帖有点粗糙，不好意思

作者: fcode 时间: 2014-4-7 10:45
你需要把 real,dimension(7)::z2,v2,q2,p2 也修改成 8 长度，否则 q2 会越界。

我不是很了解你的算法，也不知道为什么你要把他们定义为长度 7 ？

我做了如下修改：
1.把上述数组的定义从 7 改为 8
2.在 1.dat 文件最后加了一个数字 0

然后得到如下结果：

               38.046001434
               38.344001770
               39.140003204
               39.917003632
               40.278003693
               40.415004730
               40.432003021
               0.000000000

我想，除了最后一个数字，应该就正常了。

你得到怎样的结果？期望中是什么结果？

作者: 武清河 时间: 2014-4-7 12:19
本帖最后由武清河于 2014-4-7 13:40 编辑

fcode 发表于 2014-4-7 10:45
你需要把 real,dimension(7)::z2,v2,q2,p2 也修改成 8 长度，否则 q2 会越界。

我不是很了解你的算法，也 ...

我又照你的提示改了下程序：z2,v2,q2,p2 数组长度还是7，i循环从1到7，输入文件不变，运行结果基本没问题。其实这就是一个水文学上水库调洪计算的试算法，谢谢你的指拨！另外，我还有一点不懂，就是我在if块中z2(i)=z2_jiashe这条语句下面加上 p2(i)=f(z2(i))，其他不变，在输出p2数组时不能得到正确结果，然而我直接输出f(z2(i))就能得到正确结果，你能帮我解释下这是怎么回事不？修改后的代码

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

program main

real,dimension(12)::z1,v1,q1,p1

real::t=10800,z2_jiashe,q_pj,f,h,g,z2_suan

real,dimension(11)::z2,v2,q2,p2

open(1,file='1.dat')

open(2,file='2.out')

read(1,*)q1,q2

!READ*,Q1,Q2

z1(1)=38.0

p1(1)=173.9

do i=1,11

  v1(i)=g(z1(i))

  q_pj=(q1(i)+q2(i))/2

    do j=-10000,10000

      z2_jiashe=z1(i)+0.0001*j

      p2(i)=f(z2_jiashe)

      v2(i)=q_pj*t-(p1(i)+p2(i))/2*t+v1(i)

      z2_suan=h(v2(i))

         if(abs(z2_jiashe-z2_suan)<1e-4)then

            z2(i)=z2_jiashe

            v2(i)=g(z2(i))

            p2(i)=f(z2(i))

            z1(i+1)=z2(i)

            p1(i+1)=p2(i)

         endif

    enddo

enddo

do i=1,11

   write(2, '(3f20.9)'),z2(i),p2(i),v2(i)

enddo

end program main

function f(x) result(f_re)

real::x,f_re

f_re=32.93*x**2-2344*x+41710

end function f

function g(x) result(g_re)

real::x,g_re

g_re=1.911e-4*x**7.297

end function g

function h(x) result(h_re)

real::x,h_re

h_re=3.24*x**0.1369

end function h

如下

1.dat

123 Bytes, 下载次数: 1

输入文件

作者: fcode 时间: 2014-4-7 21:26
下次请直接回复，如果修改原来的帖子，别人不容易发现

我建议你学习一下单步 Debug 跟踪各变量，这样你可以自行解决大多数错误。

http://www.fcode.cn/guide-44-1.html

作者: 武清河 时间: 2014-4-8 11:53

fcode 发表于 2014-4-7 21:26
下次请直接回复，如果修改原来的帖子，别人不容易发现

我建议你学习一下单步 Debug 跟踪各变量，这样你 ...

我看过这个的，可还是找不出问题来

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)