Ada：将浮点数转换为十进制

问：

这篇文章链接到这个 Ada 2005 访问类型。目标是使用 Ada 十进制类型来获得与手动（和计算器）计算类似的结果，其中在每个中间步骤中使用了 6 位小数。

从下表可以看出，当使用欧拉方法进行进一步迭代时，使用 Ada 码获得的值开始与最后一位数字的手动计算不同。

Ada 代码的问题之一是主代码 diff.adb 中的行：return 2 * Real（XY）*;如果我也将其保留为返回 2 * X * Y，也没关系。

微分方程 （O.D.E.） 正在使用基本的欧拉方法（这是一种近似方法，不太准确）进行求解。 D.E. 为 dy/dx = 2xy。初始条件为 y0（x=x0=1） = 1。解析解为 y = e^（（x^2）-1）。目标是获得 y（x=1.5）。

我们从点 （x0，y0） = （1,1） 开始。我们使用步长 h = 0.1，即 x 随着欧拉方法中的每次迭代而增加到 1.1、1.2、1.3,..等等，y（正在寻求解的变量）的对应值由欧拉算法确定，即：

y（n） = y（n-1） + h * f（x（n-1）， y（n-1））

这里 y（n-1） 当我们启动算法时，y（0） = 1。此外，x（n-1） 是我们的起始 x（0） = 1。函数 f 是上面给出的导数函数 dy/dx = 2xy。

简而言之，h * f（x（n-1）， y（n-1）） 是“两个连续 x 值之间的水平距离”乘以梯度。梯度公式为 dy/dx = delta y /delta x，它给出 delta y 或（y 的变化）为

增量 y = 增量 x * dy/dx。

在欧拉公式中，h 是 delta x，dy/dx 是梯度。所以 h * f（x（n-1）， y（n-1）） 给出 delta y，即 delta y 的值变化，即 delta y。然后，将 y 中的这种变化添加到 y 的先前值中。欧拉方法基本上是一阶泰勒近似，x 变化很小。在曲线上绘制一条梯度线，解变量 y 的下一个值位于这条切线上，连续值为 x，即 xnew = xold + h，其中 h 是步长。

接下来的表格显示了通过手动（和计算器）计算的欧拉方法对变量 y 的解值，通过我的 Ada 代码计算，最后在最后一列中显示了精确的解。

例如，通过手工和计算器，y（x=1.1） 即 y（1） 在 x = x（1） 时计算为 y（x=1.1） = y（0） + h * f（x=1，y=1） = 1 + 0.1 * （2 * 1* 1） = 1.200000 到 6 d.p.

y（2） 在 x = x（2） 时计算为 y（x=1.2） = y（1） + h * f（x=1.1，y=1.200000） = 1.200000 + 0.1 * （2 * 1.1* 1.200000） = 1.464000 至 6 d.p.

y（3） 在 x = x（3） 时计算为 y（x=1.3） = y（2） + h * f（x=1.2，y=1.464000） = 1.464000 + 0.1 * （2 * 1.2* 1.464000） = 1.815360 至 6 d.p.

y（4） 在 x = x（4） 处计算为 y（x=1.4） = y（3） + h * f（x=1.3，y=1.815360） = 1.815360 + 0.1 * （2 * 1.3* 1.815360） = 2.287354 至 6 d.p.

y（5） 在 x = x（5） 时计算为 y（x=1.5） = y（4） + h * f（x=1.4，y=2.287354） = 2.287354 + 0.1 * （2 * 1.4* 2.287354） = 2.927813 至 6 d.p.

现在我想修改代码，以便它们使用固定数量的小数位，即小数点后 6 位。

主代码为 diff.adb：

with Ada.Text_IO;
with Euler;
procedure Diff is

  type Real is delta 0.000001 digits 9;
  type Vector is array(Integer range <>) of Real;
  type Ptr is access function (X: Real; Y: Real) return Real;
   
  package Real_IO is new Ada.Text_IO.Decimal_IO(Num => Real); 
  use Real_IO; 

  procedure Solve is new Euler(Decimal_Type => Real, Vector_Type => Vector, Function_Ptr => Ptr);

  function Maths_Func(X: Real; Y: Real) return Real is
  begin
    return 2 * Real(X*Y);
  end Maths_Func;

  
  Answer: Vector(1..6);

begin
  Solve(F => Maths_Func'Access, Initial_Value => 1.0, Increment => 0.1, Result => Answer);
  for N in Answer'Range loop
    Put(1.0 + 0.1 * Real(N-1), Exp => 0);
    Put( Answer(N), Exp => 0);
    Ada.Text_IO.New_Line;
  end loop;
end Diff;

然后是 euler.ads：

generic
  type Decimal_Type is delta <> digits <>;
  type Vector_Type is array(Integer range <>) of Decimal_Type;
  type Function_Ptr is access function (X: Decimal_Type; Y: Decimal_Type) return Decimal_Type;
procedure Euler(
  F: in Function_Ptr; Initial_Value, Increment: in Decimal_Type; Result: out Vector_Type);

和包体 euler.adb

procedure Euler
  (F : in Function_Ptr; Initial_Value, Increment : in Decimal_Type; Result : out Vector_Type)
is
   Step : constant Decimal_Type := Increment;
   Current_X : Decimal_Type := 1.0;

begin
   Result (Result'First) := Initial_Value;
   for N in Result'First + 1 .. Result'Last loop
      Result (N) := Result (N - 1) + Step * F(Current_X, Result (N - 1));
      Current_X := Current_X + Step;
   end loop;
end Euler;

在编译时，我收到指向 diff.adb 的消息：

无法从上下文中确定类型

需要显式转换为结果类型

对于行，返回 2.0 倍 X 倍 Y;

也许 2.0 在这里造成了麻烦。如何将这个浮点数转换为十进制？

我相信在diff.adb的下方，我会遇到同样的问题：

Solve(F => Maths_Func'Access, Initial_Value => 1.0, Increment => 0.1, Result => Answer);

因为它也包含浮点数。

编译是在 Windows 上使用 2011 年的 32 位 GNAT 社区版完成的。为什么是2011年？这是因为我更喜欢那一年的 IDE，而不是最近几年出现的苍白的 IDE。

接下来给出基于垃圾神代码的修订代码：

主文件 diff.adb

with Ada.Numerics.Generic_Elementary_Functions; use Ada.Numerics;
with Ada.Text_IO;                               use Ada.Text_IO;
with Euler;

procedure Diff is

   type Real is digits 7;
   type Vector is array (Positive range <>) of Real;
   type Ptr is access function (X : Real; Y : Real) return Real;
   type Round_Ptr is access function (V : Real) return Real;

   procedure Solve is new Euler (Float_Type => Real, Vector => Vector, Function_Ptr => Ptr, Function_Round_Ptr => Round_Ptr);
   package Real_Functions is new Generic_Elementary_Functions (Real);
   use Real_Functions;
   package Real_IO is new Ada.Text_IO.Float_IO (Real);
   use Real_IO;

   function DFDX (X, Y : Real) return Real is (2.0 * X * Y);
   function F (X : Real) return Real is (Exp (X**2.0 - 1.0));
   function Round (V : in Real) return Real is (Real'Rounding (1.0E6 * V) / 1.0E6);

   XI      : constant Real := 1.0;
   YI      : constant Real := 1.0;
   Step    : constant Real := 0.1;
   Result  : Vector (Positive'First .. 6); --11 if step = 0.05
   X_Value : Real;

begin
   Solve (DFDX'Access, Round'Access, XI, YI, Step, Result);
   Put_line("        x      calc     exact     delta");
   for N in Result'Range loop
      X_Value := 1.0 + Step * Real (N - 1);
      Put (X_Value, Exp => 0);
      Put (" ");
      Put (Result (N), Exp => 0);
      Put (" ");
      Put (F (X_Value), Exp => 0);
      Put (" ");
      Put (Result (N) - F (X_Value), Exp => 0);
      Ada.Text_IO.New_Line;
   end loop;
end Diff;

文件 euler.ads

generic
   type Float_Type is digits <>;
   type Vector is array (Positive range <>) of Float_Type;
   type Function_Ptr is access function (X, Y : Float_Type) return Float_Type;
   type Function_Round_Ptr is access function (V : Float_Type) return Float_Type;
procedure Euler
  (DFDX : in Function_Ptr; Round : Function_Round_Ptr; XI, YI, Step : in Float_Type; Result : out Vector);

文件 euler.adb

procedure Euler
  (DFDX : in Function_Ptr; Round : Function_Round_Ptr; XI, YI, Step : in Float_Type; Result : out Vector)
is
   H : constant Float_Type := Step;
   X : Float_Type          := XI;
begin
   Result (Result'First) := YI;
   for N in Result'First + 1 .. Result'Last loop
       Result (N) :=  Round(Result (N - 1)) + Round(H * DFDX (X, Result (N - 1)));
       X          := X + Step;
   end loop;
end Euler;

给出 **Step H = 0.1 ** 的输出

calc （Ada） 结果与手动（和计算器）计算一致。