老哥们能不能看一下我的代码为什么运行速度这么慢

鸭吼呜 · 发表于 2023-2-21 17:19:37

积分区间一大结果就跑不出来，而且精度还差的很远，答案里的代码精度很高运行速度也快，不过看不懂他这个算法

鸭吼呜 · 发表于 2023-2-21 17:20:59

这个是我的
program main
implicit none
real(kind=4),parameter :: x=1000000000
real(kind=4) i,pi,F,y,s
F=0
i=0
s=1_4/x
do while(i<x)
      y=1_4/(1_4+i*i)
      F=F+y*s
      i=i+s
end do
pi=4_4*F
write(*,*) pi
end program

鸭吼呜 · 发表于 2023-2-21 17:22:29

这个是答案里的
program pi_approx
implicit none
integer, parameter :: lk = selected_int_kind(12)
integer, parameter :: dk = selected_real_kind(13,100)
integer(lk), parameter :: nsteps = 1000000000_lk
integer(lk) :: i
real(dk) :: pi, step, sum, x

step = 1.0_dk / real(nsteps,dk)
sum = 0.0_dk
do i = 1, nsteps
x = (i - 0.5_dk) * step
sum = sum + 1.0_dk / (1.0_dk + x*x)
enddo
write(*,*) sum
pi = 4.0_dk * step * sum
write(*, '(a,f20.17)') 'approximation of pi is ', pi
end program

Transpose · 发表于 2023-2-21 18:02:37

do while(i<x)
这里应该是i<1.0，如果 i<x 计算量翻了10亿倍

你这里用的单精度，可能在很大的x的时候会损失精度(i=i+s这一步，当i和s的比值大于10**7的时候，一般加法就会损失精度)
另外他这里用的是梯形法，你这算是矩形法，按理说精度差一点

鸭吼呜 · 发表于 2023-2-21 18:44:27

Transpose 发表于 2023-2-21 18:02
do while(i

他这里好像也不是梯形法，应该是被积函数中值，但是比较费解的是，原题积分上下限是0-无穷，可是这里答案x乘了一个step，他积分区间最后应该变成了0-1，还是求出了答案，我就理解不了了，恳请大佬指点

Transpose · 发表于 2023-2-21 19:01:33

本帖最后由 Transpose 于 2023-2-21 19:03 编辑

取中间值就是梯形法

原函数是 atan(x)，x=1就是pi/4,x=无穷就是pi/2
另外这个积分区间[0,1)积分值是pi/4,最后乘以4就是pi了，没问题
如果积分区间是[0,无穷]，那么积分值是pi/2,最后就需要乘以2，但是用数值的办法比较难精确算出来

你可以做一个换元，在[1,无穷]的时候令 t=1/x ,换完之后也就变回[0,1]了

鸭吼呜 · 发表于 2023-2-21 20:10:55

Transpose 发表于 2023-2-21 19:01
取中间值就是梯形法

原函数是 atan(x)，x=1就是pi/4,x=无穷就是pi/2

哦哦对对对，我一直以为是算法上的问题没往积分上想，感激不尽

[求助] 老哥们能不能看一下我的代码为什么运行速度这么慢

规矩勋章