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

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的newton方法.但我不知道如何解决积分的上界.

根据微积分的基本定理,我认为积分的导数是积分内的函数,这是我试图重新创建的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

观察牛顿-拉夫森方法确定函数的根。即，您需要有一个函数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;