Root finding algorithm

Zero-point or root finding

PROTIDE uses a root finding algorithm to investigate each kilometer position. The root finding algorithm finds the times that the chance of touching the bottom is exactly equal to the safety criterion (this can be measured in P(bottom touch) or just a minimum UKC criterion). These roots indicate the times where safe shifts to not safe or the other way around. Therefore between these roots, on the "safe" side, we find time intervals that it is safe to be at that kilometer position.

Continuous function

A root finding algorithm needs a function it can evaluate to check if the value is positive of negative for a certain time. This function (must be continuous) is: The search method searches for the points where the above function evaluation changes sign. The method is called for a fixed location, so the only parameter is "time".

Fixed intervals and bisection

1. First with a fixed step size (15 minutes) the function is evaluated between the estimated time of arrival of the ship and 24 hours later. These 97 function evaluations and their results determine the discriminating power (or the granularity) of the algorithm. Now we are interested in the function evaluations (times) where the result changed sign. Because (under the assumption of a continuous function) we know from the intermediate value theorem that there must be a root (time where the chance of touching the bottom is exactly equal to the safety criterion) between the two points.
2. Secondly we apply a bisection algorithm on the intervals where the sign flipped to find the root with more precision. The precision is set to 2 minutes, which makes the maximum error equal to 2 minutes and the average error is estimated to 1 minute.

Bisectie

The essence of the bisection algorithm is dividing the interval in half and re-evaluating the function. This way we double the precision of our zero-point each function evaluation. The code below shows how the root is found to a certain precision in the interval:

[StartInterval.Parameter, EndInterval.Parameter]

of which we know that:

StartInterval.Result * EndInterval.Result < 0

while EndInterval.Parameter - StartInterval.Parameter > Precision / 2 do
begin
  Middle := (EndInterval.Parameter + StartInterval.Parameter) / 2;
  MiddleFunctionValue := ContinuousFunction.Execute(Middle);

  if StartInterval.Result * MiddleFunctionValue > 0 then
  begin
    StartInterval.Parameter := Middle;
    StartInterval.Result    := MiddleFunctionValue;
  end
  else
  begin
    EndInterval.Parameter := Middle;
    EndInterval.Result    := MiddleFunctionValue;    
  end;
end;

Result := (StartInterval.Parameter + EndInterval.Parameter) / 2;