# Using Ada's access type and generics to solve an ODE

I have produced Ada codes which uses the Euler algorithm Euler method to solve an ordinary differential equation. The codes given below are working fine. These codes were initially based on Ben Ari's Ada book: Ada for Software Engineers 2nd ed. 2009 Edition Section 13.6 in which Ada codes are given using the access type and generics to solve an ordinary differential equation with the Euler method. I've modified the codes so as to mimic hand (and calculator) computations so that what I get on paper matches the results I get with the Ada codes. An accuracy of six decimal places after the decimal was used.

The codes comprise of three parts: the main file diff.adb, the Euler procedure specification euler.ads and its body euler.adb. Sample results are given afterwards. Ada is robust in defining data in that which data type is more appropriate and since it's my first time with the access type and generics, I would like to know how the codes can be made better for example with the naming conventions and the data types used.

with Ada.Numerics.Generic_Elementary_Functions; use Ada.Numerics;
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);
end loop;
end Diff;



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);



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;



The output with step h = 0.1 is

1.1 1.200000 1.233678 1.233678
1.2 1.464000 1.552707 -0.033678
1.3 1.815360 1.993716 -0.088707
1.4 2.287354 2.611696 -0.178356
1.5 2.927813 3.490343 -0.562530

The calc (Ada) results agree with hand (and calculator) computations (not shown here).

Several elements merit comment:

Precision: Numerical methods using floating point math generally benefit from higher precision. There's no reason to round results during iteration. Instead, use the available precision, and format the displayed result as warranted:

type Real is digits System.Max_Digits;
Aft : constant := 6;
Put (F (X_Value), Aft => Aft, Exp => 0);


While this is a good general approach, a particular requested decimal precision may be able to rely on the Default_Aft, which is one less than the specified precision.

Initialization: For consistency, consider initializing YI with F(XI) and X_Value with XI.

Array indexing: While Ada conveniently permits indexing by any discrete type, the iterative nature of the algorithm may be more clear with a zero-based approach:

type Vector is array (Natural range <>) of Real;


Iteration: As the result Vector returns only ordinates, use the same iteration scheme to create and examine the results. Alternatively; consider passing two instances of Real_Vector or a Real_Matrix.

The compete example shown here produces the following result.

        x      calc     exact     delta
1.000000  1.000000  1.000000  0.000000
1.100000  1.200000  1.233678  0.033678
1.200000  1.464000  1.552707  0.088707
1.300000  1.815360  1.993716  0.178356
1.400000  2.287354  2.611696  0.324343
1.500000  2.927813  3.490343  0.562530
1.600000  3.806156  4.758821  0.952665
1.700000  5.024126  6.619369  1.595242
1.800000  6.732329  9.393331  2.661002
1.900000  9.155968 13.599051  4.443083
2.000000 12.635236 20.085537  7.450301

• Thank you for your points. Regarding your first one on Precision, note that the idea is not just to use Ada programming to solve an ODE. The idea is to mimic hand and calculator calculations as mentioned in the question. Books on numerical methods in their explanation will show how the solution Y is being calculated at each step and all answers are done at a fixed number of decimal places after the decimal. This is the way to explain students how the algorithm works. At exams, students are also asked to solve an ODE for instance with the Euler method on paper with say 6 d.p. precision. Mar 25, 2023 at 7:10
• Setting Aft explicitly accomplishes this, while preserving precision during iteration. Mar 25, 2023 at 12:43
• AFT I think is for display purpose only. The goal is not to preserve precision during iteration. If I am doing it on paper, Y1 may be 1.234567123....and I will only write Y1=1.234567 with 6 d.p. after the decimal point. I don't care what comes afterwards as on paper I am not preserving precision and because the question in the book say or the teacher is asking me to only keep 6 d.p. in intermediate computations. I will continue on paper with Y1 = 1.234567 in the computation of Y2 without caring that Y1 may be 1.234567123....with more precision figures. And yes I'll lose precision. Mar 25, 2023 at 13:33