Fortran双精度精度异常：输出epsilon值与原始值偏差问题

从你提供的代码和输出结果来看，epsilon值与设定值偏差较大的核心原因主要集中在变量类型不匹配和输出格式精度不足两个方面，以下是具体分析和修复方案：

一、核心问题分析 #

1. 隐式类型导致的精度损失 #

Fortran默认会根据变量首字母判断类型（首字母I-N为整数，其余为单精度实数），你的eps1变量以"e"开头，默认被声明为单精度实数（real），但你赋值的是双精度常量（1d-1、1d-3等）。单精度实数的精度只有约6-7位有效数字，无法精确存储双精度常量的数值，导致eps1中的值从一开始就存在近似误差，输出自然会偏离设定值。

比如你设定的1d-1（0.1，双精度精确值），赋值给单精度变量后会变成近似值0.099999994，用E12.2格式输出时就会被四舍五入为0.99E-01，和预期的0.10E+00不符。

2. 输出格式符适配性不足 #

你使用的E12.2格式是针对单精度实数的输出格式，仅保留两位小数的尾数部分，对于双精度变量来说，这种格式的精度不足以完整展示数值，进一步放大了显示偏差。

3. 潜在风险：参数传递与未初始化变量 #

• 如果solve_ball函数的eps形参没有声明intent(in)，默认的引用传递可能导致函数内部意外修改eps1(i)的值；
• 代码中y_max(1)和z_max(1)未初始化就参与计算，会引入未定义的数值，虽然这不是epsilon偏差的直接原因，但会影响delta-y/delta-z的结果准确性。

二、修复方案 #

1. 添加implicit none并显式声明变量类型 #

这是Fortran编程的最佳实践，彻底避免隐式类型带来的问题。在代码开头添加：

implicit none
! 显式声明所有变量的类型
double precision :: v0, t_max, h, phi, eps1(4), y_max(2), z_max(2)
integer :: i
! 如果x_max未声明，也要补充：
! double precision :: x_max

这样eps1被明确为双精度数组，赋值双精度常量时不会损失精度。

2. 修改输出格式符为双精度适配格式 #

将输出语句中的E12.2替换为ES12.4（科学计数法，保留4位有效小数）或D12.4（双精度专属格式），确保能准确显示双精度数值：

write(3,"(1i3,1i7,3ES12.4,1i6)") i, 2, eps1(i), abs(y_max(2)-y_max(1)), abs(z_max(2)-z_max(1)), ubound(r,2)

ES格式会将数值显示为1.0000E-01这样的规范科学计数法，比E格式更直观，也能保留足够精度。

3. 确保solve_ball函数的参数安全性 #

修改solve_ball函数的形参声明，添加intent(in)标记，禁止函数内部修改输入的epsilon值：

function solve_ball(h, t_max, eps, v0, phi, method) result(r)
    double precision, intent(in) :: h, t_max, eps, v0, phi
    character(len=*), intent(in) :: method
    ! 其余函数定义
end function solve_ball

4. 初始化基准值y_max(1)和z_max(1) #

在循环前，先用一个可靠的基准方法（比如更高阶的积分方法或更小的步长）计算并初始化y_max(1)和z_max(1)，确保差值计算的合理性：

! 示例：用rk4或极小步长计算基准值
r = solve_ball(h, t_max, 1d-10, v0, phi,'rk4')
y_max(1) = r(2,ubound(r,2))
z_max(1) = r(3,ubound(r,2))
deallocate(r)

三、验证效果 #

修复后，输出的epsilon值应该会和你设定的1d-1、1d-3等完全一致（或在双精度精度范围内无偏差），比如第一个epsilon会显示为1.0000E-01，而不是0.99E-01。

内容的提问来源于stack exchange，提问作者אבנר יעקב