如何在 Fortran 中获得泰勒级数的正弦
How to obtain a taylor series of sine in Fortran
我的目标是找到正弦的麦克劳林级数，其中程序 return 3 件事：用户提供的以度为单位的角度扩展，精度为 10^⁻15，第一项丢弃以获得该精度和从函数 sine(x) 计算的正弦。这是我的程序：
program serie de taylor seno
implicit none
integer*8 N, i, j
real*8 serie, x, xrad
real*8 xtest, pi, senx, func, ten, fat
write(*,*)"Choose en angle in degrees"
read(5,*) x
fat = 1.d0
ten = 10**(-15)
senx = 0.d0
pi= datan2(0.d0,-1.d0)
xrad = (x*pi)/180.d0
xtest = (x*pi)/180.d0
func =dsin(xtest)
i = 0
do while (serie .ge. ten)
do j = 1, 2*i +1
fat = fat*j
end do
N = (2*i)+1
serie = (((-1)**i)*(xrad**N))/fat
senx = senx + serie
i = i + 1
fat = 1
end do
write(*,*)senx
write (*,*) i
write(*,*)func
end program serie de taylor seno
我的问题是 do while 的循环对于任何角度只运行两次。我做错了什么？
总的来说，你的想法是正确的，但有几个问题。当尝试像这样调试时，您应该在每次循环迭代期间包括详细的输出，以尝试准确跟踪正在发生的事情。我发现的问题是：
- 您必须将公差与
serie 的 幅度 进行比较。泰勒级数的罪恶包括负项，第一个负项导致你的循环退出（在第二个，每次）。
- 您设置的公差实际上是0。
10**(-15)只包含整数，所以答案被评估为整数。使用 10.**(15)，或者更好的 1.e-15.
- 你的程序名不能有空格，至少我的
gfortran和ifort版本是这样。编辑：显然这是因为我将您的代码复制为自由形式，谢谢 francescalus
- 为了确保第一次进入
while 循环，在评估任何系列之前，将 serie 定义为某个较大的值。
此外，我还有以下其他建议：
- 不要使用不可移植的种类修饰符，例如
real*8。根据系统的不同，您会得到不同的答案。 Select 你的种类，例如，REAL_SELECTED_KIND 或 REAL64。因为您要求的是给定的宽容度，所以最好的办法是要求您的真实宽容度。
- 同样，不要使用
datan 和 dsin - 使用泛型函数将是类型无关的，并且会计算给定变量的精度。
- 不需要在循环结束时将阶乘设置为
1，这需要单独的初始化，只需在循环之前将其设置为计算阶乘即可。
- 此处不需要 8 字节整数。
- 仅仅因为最后一项小于公差并不一定意味着级数收敛到该公差。许多项的总和可能会给您带来问题，尽管对于快速收敛的 sin 级数来说这可能不是问题。
- Unit
5 通常是标准输入，但并非总是如此。 read(*,... 将连接到标准输入，或使用 ISO_FORTRAN_ENV. 中的 INPUT_UNIT
xteste有什么用？应该是 xtest?- 将不应更改的变量（例如 pi 和 tolerance）设置为参数。
- 评论你的代码！现在可读性不是很好。想一想如果你一年后再回来，你会怎么想。我还保留了用于检查代码是否正常工作的调试输出。这种输出方式适合新手学习
总而言之，您的固定程序如下所示：
program taylor
implicit none
integer, parameter :: wp = selected_real_kind(15,300)
real(wp), parameter :: pi = atan2(0._wp,-1._wp)
real(wp), parameter :: ten = 1.e-15_wp
integer :: N, i, j
real(wp) :: serie, x, xrad
real(wp) :: xtest, senx, func, fat
! -- Initialization
write(*,*) "Choose an angle in degrees"
read(*,*) x
senx = 0._wp
xrad = (x*pi)/180._wp
xtest = (x*pi)/180_wp
func = sin(xtest)
! -- Main loop
i = 0
serie = huge(serie)
do while (abs(serie) .ge. ten)
! -- Compute factorial
fat = 1._wp
do j=1,2*i+1
fat = fat*j
enddo
! -- Evaluate i-th term in series
N = (2*i)+1
serie = (((-1)**i)*(xrad**N))/fat
senx = senx + serie
i = i + 1
write(*,*) 'Completed iteration ', i
write(*,*) 'serie, senx: ', serie, senx
enddo
write(*,*) 'Completed while loop'
write(*,*) 'senx: ', senx
write (*,*) 'i: ', i
write(*,*) 'func: ', func
end program taylor
这为我提供了正确的输出：
mach5% gfortran main.f90 && ./a.out
Choose an angle in degrees
30
Completed iteration 1
serie, senx: 0.52359877559829882 0.52359877559829882
Completed iteration 2
serie, senx: -2.3924596203935038E-002 0.49967417939436376
Completed iteration 3
serie, senx: 3.2795319442867078E-004 0.50000213258879245
Completed iteration 4
serie, senx: -2.1407197692357951E-006 0.49999999186902322
Completed iteration 5
serie, senx: 8.1512566573875745E-009 0.50000000002027989
Completed iteration 6
serie, senx: -2.0315575399030637E-011 0.49999999999996431
Completed iteration 7
serie, senx: 3.5702758612702185E-014 0.50000000000000000
Completed iteration 8
serie, senx: -4.6610066605152958E-017 0.49999999999999994
Completed while loop
senx: 0.49999999999999994
i: 8
func: 0.49999999999999994
关于将整数值分配给实数的最后说明，例如 fat=1._wp。简单地说 fat=1 或 fat=1. 是一样的。但是，使用相同的表示法（例如 fat=1.5）分配非整数值会导致精度损失，因此正确的表示法是 fat=1.5_wp。所以我倾向于保持我的符号在任何地方都一致并使用 fat=1._wp 为简单起见。
我的目标是找到正弦的麦克劳林级数，其中程序 return 3 件事：用户提供的以度为单位的角度扩展，精度为 10^⁻15，第一项丢弃以获得该精度和从函数 sine(x) 计算的正弦。这是我的程序：
program serie de taylor seno
implicit none
integer*8 N, i, j
real*8 serie, x, xrad
real*8 xtest, pi, senx, func, ten, fat
write(*,*)"Choose en angle in degrees"
read(5,*) x
fat = 1.d0
ten = 10**(-15)
senx = 0.d0
pi= datan2(0.d0,-1.d0)
xrad = (x*pi)/180.d0
xtest = (x*pi)/180.d0
func =dsin(xtest)
i = 0
do while (serie .ge. ten)
do j = 1, 2*i +1
fat = fat*j
end do
N = (2*i)+1
serie = (((-1)**i)*(xrad**N))/fat
senx = senx + serie
i = i + 1
fat = 1
end do
write(*,*)senx
write (*,*) i
write(*,*)func
end program serie de taylor seno
我的问题是 do while 的循环对于任何角度只运行两次。我做错了什么？
总的来说，你的想法是正确的，但有几个问题。当尝试像这样调试时，您应该在每次循环迭代期间包括详细的输出，以尝试准确跟踪正在发生的事情。我发现的问题是：
- 您必须将公差与
serie的 幅度 进行比较。泰勒级数的罪恶包括负项，第一个负项导致你的循环退出（在第二个，每次）。 - 您设置的公差实际上是0。
10**(-15)只包含整数，所以答案被评估为整数。使用10.**(15)，或者更好的1.e-15. - 你的程序名不能有空格，至少我的
gfortran和ifort版本是这样。编辑：显然这是因为我将您的代码复制为自由形式，谢谢 francescalus - 为了确保第一次进入
while循环，在评估任何系列之前，将serie定义为某个较大的值。
此外，我还有以下其他建议：
- 不要使用不可移植的种类修饰符，例如
real*8。根据系统的不同，您会得到不同的答案。 Select 你的种类，例如，REAL_SELECTED_KIND或REAL64。因为您要求的是给定的宽容度，所以最好的办法是要求您的真实宽容度。 - 同样，不要使用
datan和dsin- 使用泛型函数将是类型无关的，并且会计算给定变量的精度。 - 不需要在循环结束时将阶乘设置为
1，这需要单独的初始化，只需在循环之前将其设置为计算阶乘即可。 - 此处不需要 8 字节整数。
- 仅仅因为最后一项小于公差并不一定意味着级数收敛到该公差。许多项的总和可能会给您带来问题，尽管对于快速收敛的 sin 级数来说这可能不是问题。
- Unit
5通常是标准输入，但并非总是如此。read(*,...将连接到标准输入，或使用ISO_FORTRAN_ENV. 中的 xteste有什么用？应该是xtest?- 将不应更改的变量（例如 pi 和 tolerance）设置为参数。
- 评论你的代码！现在可读性不是很好。想一想如果你一年后再回来，你会怎么想。我还保留了用于检查代码是否正常工作的调试输出。这种输出方式适合新手学习
INPUT_UNIT
总而言之，您的固定程序如下所示：
program taylor
implicit none
integer, parameter :: wp = selected_real_kind(15,300)
real(wp), parameter :: pi = atan2(0._wp,-1._wp)
real(wp), parameter :: ten = 1.e-15_wp
integer :: N, i, j
real(wp) :: serie, x, xrad
real(wp) :: xtest, senx, func, fat
! -- Initialization
write(*,*) "Choose an angle in degrees"
read(*,*) x
senx = 0._wp
xrad = (x*pi)/180._wp
xtest = (x*pi)/180_wp
func = sin(xtest)
! -- Main loop
i = 0
serie = huge(serie)
do while (abs(serie) .ge. ten)
! -- Compute factorial
fat = 1._wp
do j=1,2*i+1
fat = fat*j
enddo
! -- Evaluate i-th term in series
N = (2*i)+1
serie = (((-1)**i)*(xrad**N))/fat
senx = senx + serie
i = i + 1
write(*,*) 'Completed iteration ', i
write(*,*) 'serie, senx: ', serie, senx
enddo
write(*,*) 'Completed while loop'
write(*,*) 'senx: ', senx
write (*,*) 'i: ', i
write(*,*) 'func: ', func
end program taylor
这为我提供了正确的输出：
mach5% gfortran main.f90 && ./a.out
Choose an angle in degrees
30
Completed iteration 1
serie, senx: 0.52359877559829882 0.52359877559829882
Completed iteration 2
serie, senx: -2.3924596203935038E-002 0.49967417939436376
Completed iteration 3
serie, senx: 3.2795319442867078E-004 0.50000213258879245
Completed iteration 4
serie, senx: -2.1407197692357951E-006 0.49999999186902322
Completed iteration 5
serie, senx: 8.1512566573875745E-009 0.50000000002027989
Completed iteration 6
serie, senx: -2.0315575399030637E-011 0.49999999999996431
Completed iteration 7
serie, senx: 3.5702758612702185E-014 0.50000000000000000
Completed iteration 8
serie, senx: -4.6610066605152958E-017 0.49999999999999994
Completed while loop
senx: 0.49999999999999994
i: 8
func: 0.49999999999999994
关于将整数值分配给实数的最后说明，例如 fat=1._wp。简单地说 fat=1 或 fat=1. 是一样的。但是，使用相同的表示法（例如 fat=1.5）分配非整数值会导致精度损失，因此正确的表示法是 fat=1.5_wp。所以我倾向于保持我的符号在任何地方都一致并使用 fat=1._wp 为简单起见。