基于Ada 2005访问类型的欧拉法求解微分方程代码问题求助

我尝试用《Ada for Software Engineers》（Ben Ari第二版）第13.6.1节（第263页起）介绍的访问类型技术，通过欧拉法求解微分方程 dy/dx = 2*x*y，初始条件为 y(x=1)=1，步长 h=0.1，需要输出 y(1) 到 y(1.5) 的值。编写了3个Ada文件：主文件diff.adb、euler.ads和euler.adb。编译时曾出现错误：euler.adb:6:36 missing argument for parameter Y，参考建议调整代码后编译通过，但运行结果和手动计算的预期值不符，现寻求代码修正方案。

现有代码 #

diff.adb

with Ada.Text_IO;
with Euler;
procedure Diff is

  type Real is digits 6;
  type Vector is array(Integer range <>) of Real;
  type Ptr is access function (X: Real; Y: Real) return Real;

  procedure Solve is new Euler(Real, Vector, Ptr);

  function Ident(X: Real; Y: Real) return Real is
  begin
    return 2.0*X*Y;
  end Ident;

  package Real_IO is new Ada.Text_IO.Float_IO(Real);
  use Real_IO;

  Answer: Vector(1..6);

begin
  Solve(Ident'Access, 1.0, 0.1, 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

--
-- Solving a differential equation.
-- Demonstrates generic floating point type.
--
generic
  type Float_Type is digits <>;
  type Vector is array(Integer range <>) of Float_Type;
  type Function_Ptr is access function (X: Float_Type; Y: Float_Type) return Float_Type;
procedure Euler(
  F: in Function_Ptr; Init, H: in Float_Type; Result: out Vector);

修改后的euler.adb

procedure Euler
  (F : in Function_Ptr; Init, H : in Float_Type; Result : out Vector)
is
      Step : constant Float_Type := H;
     
     Current_X : Float_Type := 1.0;
begin
   Result (Result'First) := Init;
   for N in Result'First + 1 .. Result'Last loop
      Current_X := Current_X + Step;
      Result (N) := Result (N - 1) + Step * F (Current_X, Result (N - 1));
   end loop;
end Euler;

预期结果 #

实际运行结果 #

代码修正方案 #

问题出在欧拉法的迭代逻辑上。标准欧拉法的公式是 y(n+1) = y(n) + h * f(x(n), y(n))，也就是用当前点的x值x(n)代入微分方程计算斜率，再更新y值。但现有代码中，先把Current_X更新到了x(n+1)，再用这个新x值去计算f，导致每一步的斜率都偏大，结果偏离预期。

修正后的euler.adb需要调整Current_X的更新顺序：先基于当前x(n)计算y(n+1)，再把x更新为x(n+1)。

修正后的euler.adb代码

procedure Euler
  (F : in Function_Ptr; Init, H : in Float_Type; Result : out Vector)
is
  Step : constant Float_Type := H;
  Current_X : Float_Type := 1.0;
begin
  Result (Result'First) := Init;
  for N in Result'First + 1 .. Result'Last loop
    -- 先基于当前x(n)计算y(n+1)
    Result (N) := Result (N - 1) + Step * F (Current_X, Result (N - 1));
    -- 再更新x到x(n+1)
    Current_X := Current_X + Step;
  end loop;
end Euler;

编译运行修正后的代码，就能得到和预期一致的结果。

内容的提问来源于stack exchange，提问作者Adaenthusiast