Fortran Coder

查看: 17014|回复: 6
打印 上一主题 下一主题

[求助] 这究竟是算法不对还是语法错了?求指点!

[复制链接]

4

帖子

1

主题

0

精华

新人

F 币
27 元
贡献
14 点
跳转到指定楼层
楼主
发表于 2014-4-6 22:24:59 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
本帖最后由 武清河 于 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

输入文件

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2033

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1641 元
贡献
709 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

沙发
发表于 2014-4-7 08:16:27 | 只看该作者
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
后面的各种越界。

4

帖子

1

主题

0

精华

新人

F 币
27 元
贡献
14 点
板凳
 楼主| 发表于 2014-4-7 10:37:43 | 只看该作者
fcode 发表于 2014-4-7 08:16
1.请给出输入文件(或其内容)。1.dat
2.建议:写代码要搞好排版,缩进。
3.目测你的代码错误很多,不知道 ...

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

2033

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1641 元
贡献
709 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

地板
发表于 2014-4-7 10:45:56 | 只看该作者
你需要把 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

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

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


4

帖子

1

主题

0

精华

新人

F 币
27 元
贡献
14 点
5#
 楼主| 发表于 2014-4-7 12:19:58 | 只看该作者
本帖最后由 武清河 于 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

输入文件

2033

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1641 元
贡献
709 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

6#
发表于 2014-4-7 21:26:58 | 只看该作者
下次请直接回复,如果修改原来的帖子,别人不容易发现

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

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

4

帖子

1

主题

0

精华

新人

F 币
27 元
贡献
14 点
7#
 楼主| 发表于 2014-4-8 11:53:45 | 只看该作者
fcode 发表于 2014-4-7 21:26
下次请直接回复,如果修改原来的帖子,别人不容易发现

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

我看过这个的,可还是找不出问题来
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-12-24 01:37

Powered by Tencent X3.4

© 2013-2024 Tencent

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