Как решить для верхнего предела интеграла, используя метод Ньютона?

Question

Я хочу написать программу, которая использует метод Ньютона:

Чтобы оценить x этого интеграла:

Где X - общее расстояние.

У меня есть функции для расчета времени, необходимого для достижения определенного расстояния, с использованием метода трапеции длячисленное интегрирование.Без использования 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

и он выбирает свои значения для скорости через ppval кубической сплайн-интерполяции между набором точек данных.Где экстраполированные значения не должны быть извлекаемыми.

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

График сплайна скорости, если вам интересно, оценивается по:

dx= 1:0.1:65

Теперь я хочу написать функцию, которая может решитьдля расстояния, пройденного через определенное время, используя метод Ньютона без fzero / fsolve.Но я понятия не имею, как решить для верхней границы интеграла.

В соответствии с фундаментальной теоремой исчисления я предполагаю, что производная интеграла является функцией внутри интеграла, что я и пыталсячтобы воссоздать как Time_to_destination / (1 / скорость), я добавил постоянную, которую я хочу найти, к времени до места назначения, поэтому ее

(Time_to_destination - (время ввода)) / (1 / скорость)

Не уверен, что я делаю это правильно.

РЕДАКТИРОВАТЬ: переписал мой код, теперь работает лучше, но мое условие остановки для Ньютона Рафсона, похоже, не сходится к нулю.Я также пытался реализовать ошибку из интеграции трапеции (ET), но не уверен, стоит ли мне ее реализовывать.Также найдите файл маршрута внизу.

Условие остановки и расчет ошибки метода Ньютона:

Оценка ошибки трапеции:

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

Файл маршрута: https://drive.google.com/open?id=18GBhlkh5ZND1Ejh0Muyt1aMyK4E2XL3C

Answer 1

Одна из основных проблем с моим кодом заключалась в том, что n было слишком низким, ошибка трапецеидальной суммы, оценка моего интеграла, была слишком высокой, чтобы метод Ньютона Рафсона мог сходиться к очень маленькому числу.

Вот мой окончательный код для этой проблемы:

function x = distance(T, route)
load(route)
n=10e6;
x = mean(distance_km);
i = 1;
maxiter=100;
tol= 5e-4;
condition=inf
fx = @(x) time_to_destination(x, route,n); 
dfx = @(x) 1./velocity(x, route);

while  condition > tol && i <= maxiter
    i = i + 1 ;
    Guess2 = x - ((fx(x) - T)/(dfx(x)))
    condition = abs(Guess2 - x)
    x = Guess2;

end
end

Answer 2

Обратите внимание, что метод Ньютона-Рафсона определяет корни функции.Т.е. у вас должна быть функция f (x) такая, что f (x) = 0 при желаемом решении.

В этом случае вы можете определить f as

f (x) = Время (x) - t

, где t - желаемое время.Тогда по второй фундаментальной теореме исчисления

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

С этими определенными функциями реализация становится довольно простой!

Сначала мы определим простую функцию Ньютона-Рафсона, которая принимает анонимные функции в качестве аргументов ( f и f '), а также первоначальное предположение 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

Затем мы можем вызвать нашу функцию следующим образом

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

Результат

>> distance
    distance = 25.5877

---

Кроме того, я переписал ваши time_to_destination иvelocity работает следующим образом.Эта версия time_to_destination использует все доступные данные для более точной оценки интеграла.При использовании этих функций метод, по-видимому, сходится быстрее.

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

Использование этих новых функций требует незначительного изменения определений анонимных функций.

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

Поскольку интеграл оценивается болееточно решение немного отличается

>> distance
    distance = 25.7771

---

Редактировать

Обновленное условие остановки может быть реализовано как незначительное изменение функции newton_method.Мы не должны ожидать, что ошибка правила трапеции обнулится, поэтому я опускаю это.

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

Чтобы проверить наш ответ, мы можем построить график зависимости времени от расстояния и убедиться, что наша оценка падает на кривой.

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