How to solve for the upper limit of an integral using Newton's method?

Question

I want to write a program that makes use of Newtons Method:

To estimate the x of this integral:

Where X is the total distance.

I have functions to calculate the Time it takes to arrive at a certain distance by using the trapezoid method for numerical integration. Without using trapz.

function T = time_to_destination(x, route, n)
h=(x-0)/n;
dx = 0:h:x;
y = (1./(velocity(dx,route)));

Xk = dx(2:end)-dx(1:end-1); 
Yk = y(2:end)+y(1:end-1);

T = 0.5*sum(Xk.*Yk);
end

and it fetches its values for velocity, through ppval of a cubic spline interpolation between a set of data points. Where extrapolated values should not be fetcheable.

function [v] = velocity(x, route)
load(route);

if all(x >= distance_km(1))==1 & all(x <= distance_km(end))==1
    estimation = spline(distance_km, speed_kmph);
    v = ppval(estimation, x);
else
    error('Bad input, please choose a new value')
end
end

Plot of the velocity spline if that's interesting to you evaluated at:

dx= 1:0.1:65

Now I want to write a function that can solve for distance travelled after a certain given time, using newton's method without fzero / fsolve . But I have no idea how to solve for the upper bound of a integral.

According to the fundamental theorem of calculus I suppose the derivative of the integral is the function inside the integral, which is what I've tried to recreate as Time_to_destination / (1/velocity) I added the constant I want to solve for to time to destination so its

(Time_to_destination - (input time)) / (1/velocity)

Not sure if I'm doing that right.

EDIT: Rewrote my code, works better now but my stopcondition for Newton Raphson doesnt seem to converge to zero. I also tried to implement the error from the trapezoid integration ( ET ) but not sure if I should bother implementing that yet. Also find the route file in the bottom.

Stop condition and error calculation of Newton's Method:

Error estimation of trapezoid:

Function x = distance(T, route)
n=180
route='test.mat'
dGuess1 = 50;
dDistance = T;
i = 1;
condition = inf;

while  condition >= 1e-4  && 300 >= i
    i = i + 1 ;
    dGuess2 = dGuess1 - (((time_to_destination(dGuess1, route,n))-dDistance)/(1/(velocity(dGuess1, route))))
    if i >= 2
        ET =(time_to_destination(dGuess1, route, n/2) - time_to_destination(dGuess1, route, n))/3;
        condition = abs(dGuess2 - dGuess1)+ abs(ET);
    end
    dGuess1 = dGuess2;
end
x = dGuess2

Route file: https://drive.google.com/open?id=18GBhlkh5ZND1Ejh0Muyt1aMyK4E2XL3C

Answer 1

Observe that the Newton-Raphson method determines the roots of the function. I.e. you need to have a function f(x) such that f(x)=0 at the desired solution.

In this case you can define f as

f(x) = Time(x) - t

where t is the desired time. Then by the second fundamental theorem of calculus

f'(x) = 1/Velocity(x)

With these functions defined the implementation becomes quite straightforward!

First, we define a simple Newton-Raphson function which takes anonymous functions as arguments (f and f') as well as an initial guess x0.

function x = newton_method(f, df, x0)
    MAX_ITER = 100;
    EPSILON = 1e-5;

    x = x0;
    fx = f(x);
    iter = 0;
    while abs(fx) > EPSILON && iter <= MAX_ITER
        x = x - fx / df(x);
        fx = f(x);
        iter = iter + 1;
    end
end

Then we can invoke our function as follows

t_given = 0.3;   % e.g. we want to determine distance after 0.3 hours.
n = 180;
route = 'test.mat';
f = @(x) time_to_destination(x, route, n) - t_given; 
df = @(x) 1/velocity(x, route);

distance_guess = 50;
distance = newton_method(f, df, distance_guess);

Result

>> distance
    distance = 25.5877

---

Also, I rewrote your time_to_destination and velocity functions as follows. This version of time_to_destination uses all the available data to make a more accurate estimate of the integral. Using these functions the method seems to converge faster.

function t = time_to_destination(x, d, v)
    % x is scalar value of destination distance
    % d and v are arrays containing measured distance and velocity
    % Assumes d is strictly increasing and d(1) <= x <= d(end)
    idx = d < x;
    if ~any(idx)
        t = 0;
        return;
    end
    v1 = interp1(d, v, x);
    t = trapz([d(idx); x], 1./[v(idx); v1]);
end

function v = velocity(x, d, v)
    v = interp1(d, v, x);
end

Using these new functions requires that the definitions of the anonymous functions are changed slightly.

t_given = 0.3;   % e.g. we want to determine distance after 0.3 hours.
load('test.mat');
f = @(x) time_to_destination(x, distance_km, speed_kmph) - t_given; 
df = @(x) 1/velocity(x, distance_km, speed_kmph);

distance_guess = 50;
distance = newton_method(f, df, distance_guess);

Because the integral is estimated more accurately the solution is slightly different

>> distance
    distance = 25.7771

---

Edit

The updated stopping condition can be implemented as a slight modification to the newton_method function. We shouldn't expect the trapezoid rule error to go to zero so I omit that.

function x = newton_method(f, df, x0)
    MAX_ITER = 100;
    TOL = 1e-5;

    x = x0;
    iter = 0;
    dx = inf;
    while dx > TOL && iter <= MAX_ITER
        x_prev = x;
        x = x - f(x) / df(x);
        dx = abs(x - x_prev);
        iter = iter + 1;
    end
end

To check our answer we can plot the time vs. distance and make sure our estimate falls on the curve.

...
distance = newton_method(f, df, distance_guess);

load('test.mat');
t = zeros(size(distance_km));
for idx = 1:numel(distance_km)
    t(idx) = time_to_destination(distance_km(idx), distance_km, speed_kmph);
end
plot(t, distance_km); hold on;
plot([t(1) t(end)], [distance distance], 'r');
plot([t_given t_given], [distance_km(1) distance_km(end)], 'r');
xlabel('time');
ylabel('distance');
axis tight;