This post is linked to this one Ada 2005 access type. The goal is to use **Ada decimal type to get similar results as to hand (and calculator) computations in which 6 decimal places have been used in each intermediate step.**

As can be seen from the table below, the **values obtained with the Ada code starts to differ from the hand calculation in the last digit when further iterations with the Euler method are taken**.

One of the issues with the Ada code was with the line in the main code **diff.adb**: *return 2 * Real(X*Y)*; It doesn't matter if I leave it as *return 2 * X * Y* as well.

The differential equation (O.D.E.) is being solved using the basic Euler method (which is an approximate method which is not that accurate).
The D.E. is dy/dx = 2*x*y. The initial condition is at y0(x=x0=1) = 1. The analytical solution is y = e^((x^2)-1). The objective is to obtain y(x=1.5).

We start with the point (x0,y0) = (1,1). We use a step size h = 0.1 i.e. x is increased with each iteration in the Euler method to 1.1, 1.2, 1.3,..etc. and the corresponding value of y (the variable whose solution is being sought) is determined from the Euler algorithm which is:

y(n) = y(n-1) + h * f(x(n-1), y(n-1))

Here y(n-1) when we start the algorithm is y(0) = 1. Also x(n-1) is our starting x(0) = 1. The function f is the derivative function dy/dx given above as dy/dx = 2*x*y.

Briefly, h * f(x(n-1), y(n-1)) is the "horizontal distance between two successive x values" multiplied by the gradient. The gradient formula is dy/dx = delta y /delta x which gives delta y or (the change in y) as

delta y = delta x * dy/dx.

In the Euler formula h is the delta x and dy/dx is the gradient. So h * f(x(n-1), y(n-1)) gives delta y which is the change in the value of y i.e. delta y. This change in y is then added to the previous value of y. The Euler method is basically a first order Taylor approximation with a small change in x. A gradient line is drawn to the curve and the next value of the solution variable y is on this tangent line at the successive value of x i.e. xnew = xold + h where h is the step.

The table next shows the solution values for the variable y by the Euler method when calculated by hand (and calculator), by my Ada code and finally in the last column the exact solution.

x |
y (hand) |
Ada code |
y (exact) |

1.1 |
1.200000 |
1.200000 |
1.233678 |

1.2 |
1.464000 |
1.464000 |
1.552707 |

1.3 |
1.815360 |
1.815360 |
1.993716 |

1.4 |
2.287354 |
2.287353 |
2.611696 |

1.5 |
2.927813 |
2.927811 |
3.490343 |

By hand and calculator for instance, y(x=1.1) i.e y(1) at x = x(1) is calculated as
y(x=1.1) = y(0) + h * f(x=1,y=1) = 1 + 0.1 * (2 * 1* 1) = 1.200000 to 6 d.p.

y(2) is calculated at x = x(2) as
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 to 6 d.p.

y(3) is calculated at x = x(3) as
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 to 6 d.p.

y(4) is calculated at x = x(4) as
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 to 6 d.p.

y(5) is calculated at x = x(5) as
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 to 6 d.p.

Now I want to modify the codes so that they work with a fixed number of decimal places which is 6 here after the decimal place.

The main code is **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;
```

Then comes **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);
```

and the package body **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;
```

On compilation, I get the messages pointing to **diff.adb**:

**type cannot be determined from context**

**explicit conversion to result type required**

for the line **return 2.0 times X times Y**;

Perhaps the **2.0** is causing the trouble here. How to convert this Float number to Decimal?

I believe that further down in **diff.adb**, I will get the same issue with the line:

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

for it contains Floating point numbers as well.

The compilation was done on Windows with the 32-bit GNAT community edition of year 2011. Why 2011? This is because I like the IDE better for that year rather than the pale ones which come in the recent years.

**The revised codes based on trashgod codes which work are given next:**

The main file **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;
```

The file **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);
```

The file **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;
```

giving the output with **step h = 0.1 **

x |
calc (Ada) |
exact |
delta |

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.