如何使用牛顿法求解积分的上限？

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/velocity) 我添加了我想要解决的常数到目的地的时间所以它

(Time_to_destination - (输入时间)) / (1/速度)

不确定我这样做是否正确。

编辑：重写了我的代码，现在工作得更好了，但我对 Newton Raphson 的停止条件似乎没有收敛到零。我还尝试从梯形积分 ( 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

观察到 Newton-Raphson 方法确定了函数的根。 IE。您需要有一个函数 f(x) 使得 f(x)=0 在所需的解决方案中。

在这种情况下，您可以将 f 定义为

f(x) = 时间(x) - t

其中 t 是所需的时间。然后由微积分第二基本定理

f'(x) = 1/速度(x)

定义了这些函数后，实现变得非常简单！

首先，我们定义一个简单的 Newton-Raphson 函数，它将匿名函数作为参数（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;

Answer 2

我的代码的一个主要问题是 n 太低，梯形和的误差，我的积分估计，对于 newton raphson 方法收敛到一个非常小的数来说太高了。

这是我解决这个问题的最终代码：

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

如何使用牛顿法求解积分的上限？

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

matlab

numerical-integration