24小时热门版块排行榜

返回列表

当前主题已经存档。

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

jianchaoyv

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 964.9
帖子: 198
在线: 94.8小时
虫号: 657809
注册: 2008-11-19
专业: 通信理论与系统

[交流] [求助]一个关于梯形积分的程序

program D3R1
!driver for routine TRAPZD
real,parameter::nmax=15,pio2=1.5707963
external Func,Fint
real::A=0.0,B=pio2
write(*,'(1x,a)') 'integral of Func with 2^(n-1) points'
write(*,'(1x,a,f10.6)')'actual value of integral is',&
            Fint(B)-Fint(A)
write(*,'(1x,t7,a,t16,a)')'n','Approx.integral'
do i=1,nmax
call TRAPZD(Func,A,B,s,i)
write(*,'(1x,i6,f20.6)')i,s
end do
end program D3R1
Function Func(x)
Func=(x**2)*(x**2-2.0)*sin(x)
end function Func
function Fint(x)
!integral of Func
Fint=4.0*x*((x**2)-7.0)*sin(x)-((x**4)-14.0*(x**2)+28.0)*cos(x)
end function Fint

subroutine TRAPZD(Func,A,B,s,n)
integer::n,tnm
integer,parameter::p=selected_real_kind(6,20)
real(p)::del,sum,x
if(n==1)then
  s=0.5*(b-a)*(Func(a)+Func(b))
  else
  tnm=2**(n-1)

  del=(b-a)/tnm
  x=A
  sum=0.0
  do j=2,tnm
x=x+del
sum=sum+Func(x)

end do
s=0.5*(Func(A)+Func(B)+2*sum)*del
end if
end subroutine TRAPZD
运算的结果如下:
integral of Func with 2^(n-1) points
actual value of integral is -0.479158
   n       Approx.integral
   1          0.905772
   2          -0.020945
   3          -0.361461
   4          -0.449584
   5          -0.471756
   6          -0.477307
   7          -0.478697
   8          -0.479042
   9          -0.479126
   10          -0.479146
   11          -0.479151
   12          -0.479152
   13          -0.479153
   14          -0.479411
   15          -0.479676
Press any key to continue
请问n=14,15 时怎么越不接近真实值?

回复此楼

» 猜你喜欢

假如你的研究生提出不合理要求已经有12人回复
实验室接单子已经有7人回复
全日制（定向）博士已经有5人回复
萌生出自己或许不适合搞科研的想法，现在跑or等等看？已经有4人回复
Materials Today Chemistry审稿周期已经有4人回复
参与限项已经有3人回复
对氯苯硼酸纯化已经有3人回复
所感已经有4人回复
要不要辞职读博？已经有7人回复
北核录用已经有3人回复

1楼 2009-03-29 17:13:19

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

程序强帖: 16
应助: 157 (高中生)
贵宾: 0.02
金币: 18844.7
红花: 29
帖子: 3803
在线: 1422.4小时
虫号: 183750
注册: 2006-02-13
专业: 污染生态化学

俺不懂 f90，写了一个 f77 的程序供你参考。我怀疑是数值类型的问题……

CODE:

C234567

      SUBROUTINE TX(FUNC,A,B,N,S)

      IMPLICIT NONE

      INTEGER N, TNM, I

      DOUBLE PRECISION DEL, SUM, X, A, B, S, FUNC

      EXTERNAL FUNC

      IF( N .EQ. 1) THEN

        S = 0.5D0 * (B - A) * (FUNC(A) + FUNC(B))

      ELSE

        TNM = 2**(N - 1)

        DEL = (B - A) / TNM

        X = A

        SUM = 0.0D0

        DO 5 I = 2, TNM

          X = X + DEL

          SUM = SUM + FUNC(X)

5       CONTINUE

        S=0.5D0*(FUNC(A) + FUNC(B) + 2D0 * SUM) * DEL

      ENDIF

      RETURN

      END

      DOUBLE PRECISION FUNCTION FUNC(X)

      IMPLICIT NONE

      DOUBLE PRECISION X

      FUNC = (X**2) * (X**2 - 2.0D0) * SIN(X)

      RETURN

      END

      PROGRAM MAIN

      IMPLICIT NONE

      INTEGER N, I

      DOUBLE PRECISION S, A, B, FUNC

      EXTERNAL FUNC

      A = 0.0D0

      B = 1.5707963D0

      N = 20

      DO 5 I=1,N

        CALL TX(FUNC,A,B,I,S)

        WRITE(*,'(I4,F20.10)') I, S

5     CONTINUE

      END

运行结果：

CODE:

   1        0.9057727803

   2       -0.0209449905

   3       -0.3614613180

   4       -0.4495837621

   5       -0.4717563216

   6       -0.4773076746

   7       -0.4786960160

   8       -0.4790431327

   9       -0.4791299138

  10       -0.4791516092

  11       -0.4791570331

  12       -0.4791583890

  13       -0.4791587280

  14       -0.4791588128

  15       -0.4791588339

  16       -0.4791588392

  17       -0.4791588406

  18       -0.4791588409

  19       -0.4791588410

  20       -0.4791588410

[ Last edited by snoopyzhao on 2009-4-1 at 13:19 ]

赞一下

回复此楼

4楼2009-04-01 11:54:19

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

查看全部 6 个回答

abingchem

木虫 (著名写手)

应助: 2 (幼儿园)
金币: 4871
散金: 968
红花: 2
帖子: 1672
在线: 117.6小时
虫号: 206219
注册: 2006-03-03
性别: GG
专业: 化工系统工程

★ ★ ★
sunxiao(金币+3,VIP+0):欢迎参与，加分鼓励 3-30 22:39

俺不会Fortran,对你的程序只能看了个大概.循环部分没有问题,至少写成matlab后执行是正确的:
clear;
clc
b=pi/2
a=0
nse=[];
sse=[];

for n=10:20
tnm=2^(n-1)

del=(b-a)/tnm
x=a
s=0.0

for j=2:tnm
x=x+del;
s=s+Func(x);
end
s=0.5*(Func(a)+Func(b)+2*s)*del
nse=[nse n];
sse=[sse s];
end
[nse'  sse']
---------------------------------------------
结果是这样的:
10.00000000000000  -0.47915157829510
  11.00000000000000  -0.47915700215368
  12.00000000000000  -0.47915835811872
  13.00000000000000  -0.47915869710992
  14.00000000000000  -0.47915878185820
  15.00000000000000  -0.47915880304466
  16.00000000000000  -0.47915880834252
  17.00000000000000  -0.47915880966840
  18.00000000000000  -0.47915881000004
  19.00000000000000  -0.47915881008306

有点怀疑的是你的变量s是在哪里初始化的,这个值是应该先初始化为零吧.

赞一下(10人)

回复此楼

冰是从最寒冷的那天开始融化的