首先，这是一个可复制的解决方案。我的答案剩下的部分是解释 :)
(你需要从 FEX 下载 FindRealRoots 才能运行此程序)。
trange = [tstart tend]
a = m2^m5;
b = m3^m6;
c = m4^m7;
X = @(t) t(:).' + m8;
g = @(t) 1 + [a b c] * bsxfun(@power, X(t), [m5; m6; m7]);
gp = @(t) ([m5 m6 m7].*[a b c]) * bsxfun(@power, X(t), [m5-1; m6-1; m7-1]);
gpp = @(t) (([m5 m6 m7]-1).*[m5 m6 m7].*[a b c]) * ...
bsxfun(@power, X(t), [m5-1; m6-1; m7-1]);
sgn = @(t) m1*(2*(gp(t)).^2 - g(t).^2*gpp(t) );
y = @(t) m1*(1 + 1./g(t));
yp = @(t) -m1*gp(t) ./ g(t).^2;
R = FindRealRoots(gp, tstart , tend, 2*ceil(max([m5 m6 m7])));
y0 = y(tstart);
yend = y(tend);
if ~isempty(R)
extrema = y(R);
signs = sgn(R);
[maximum, Imax] = max(extrema(signs < 0));
tmax = R(Imax);
if isempty(tmax)
[~,ind] = max([y0, yend]);
tmax = trange(ind);
end
ymax = y(tmax);
[minima, Imin] = extrema(signs > 0);
tmin = R(Imin);
if isempty(tmin)
[~,ind] = min([y0, yend]);
tmin = trange(ind);
end
if any(tmin < tmax)
tmin = tmin(find(tmin < tmax, 1, 'last'));
else
[~,ind] = min([y0, yend]);
tmin = trange(ind);
end
ymin = y(tmin);
else
[ymax,ind] = max([y0, yend]);
tmax = trange(ind);
[ymin,ind] = min([y0, yend]);
tmin = trange(ind);
end
现在，解释一下。
虽然Dan的答案是一个不错的第一步，但我觉得我必须在一些地方进行更正。
首先，数值导数只在以下情况下按照Dan所描述的方式计算：
- 您有网格数据（在许多不同点t处的y(t)的值）
- 已知底层模型函数y(t)不允许插值
而这在您的情况下并不正确。当您有一个可以在任意点评估的显式函数时，最好使用
有限差分来计算数值导数，通常仅在无法通过分析方法找到显式导数时使用。
例如，在x=π/4处（使用h=0.001的中心差分）的cos(x)的数值导数：
z'(π/4) ≅ ( z(x+h)-z(x-h) ) / 2h
= ( cos(π/4 + 0.001) - cos(π/4 - 0.001) ) / 0.002
= -0.7071066633353995...
-sin(π/4) = -0.7071067811865475... (value of analytical derivative)
说了这些，我们转向更重要的一点；在您不得不求助于数值方法之前，可以将数学推进得比Dan更远。这样做可能会更加繁琐，但它将改善您对问题、MATLAB、数学的理解，并提高结果的准确性和稳健性。
对于您的函数
y = f(t)，您声明有三个目标：
- 找到速度曲线
- 找到最大值
y(tmax) = ymax
- 找到发生在最大值
ymax 之前的最小值 ymin，tmin
当翻译成数学术语时：
- 求导数
y' = dy/dt
- (和3.) 找出所有
y' 的根 R，并在所有 R 处计算 y 和 y''，还要在初始时间 t0 处计算 y。
让我们从目标1开始。您有一个关于时间的显式函数 y = f(t)，这基本上是有理函数的一种变形形式。
y(t) = a + f(t)/g(t)
为了简洁起见，我们可以省略（t）部分。这样写起来更容易，而且可以理解y、f和g都依赖于t。所以，
y = a + f/g
你有一个优势，即
f(t) = <constant> = a = m1，这在后面的事情中会更加方便。这种函数的导数表达式采用以下形式（即
商规则）：
derivative。
y' = (f'·g - f·g') / g²
由于 f(t) = m1 = <constant>，所以它对时间的导数为零：
y' = -m1·g' / g²
现在，我们需要找到
g的导数。函数
g是多项式：
g = 1 + (m2·(t+m8))m5+ (m3·(t+m8))m6+ (m4·(t+m8))m7
让我们定义
X = t + m8，
a = m2m5，
b = m3m6和
c = m4m7。这样，
g看起来会少吓人很多：
g = 1 + a·Xm5+ b·Xm6+ c·Xm7
Differentiation of polynomials follows the
elementary power rule (and
chain rule, but that's not
really relevant here because
X' = 1):
g' = m5·a·Xm5-1+ m6·b·Xm6-1+ m7·c·Xm7-1
因此，
y' = -m1·g' / g² = -m1 · (m5·a·Xm5-1+ m6·b·Xm6-1+ m7·c·Xm7-1) / (1 + a·Xm5+ b·Xm6+ c·Xm7)²
目标1已完成！您现在可以在任何地方使用完全准确的速度曲线。下面是如何在MATLAB中实现此函数：
a = m2^m5;
b = m3^m6;
c = m4^m7;
X = @(t) t + m8;
dydt = @(t) -m1*( [m5*a m6*b m7*c]*X(t).^[m5-1; m6-1; m7-1] ) ./ ...
( 1 + [a b c]*X(t).^[m5; m6; m7] ).^2;
现在是目标2：
- 找到
y'=0的所有根。
- 在这些根上评估
y。
- 在这些根上评估
y''。
- 在
t区间的端点上评估y。
步骤3是必要的，以确定y'=0的根R是否描述了y的最大值或最小值。如果y''(R) < 0，则处理的是最大值，如果y'' > 0，则处理的是最小值，如果y''(R) = 0，则处理的是一个拐点，可以安全地忽略此问题。
第四步是必要的，因为在函数的最大值之前，可能存在没有最小值的情况，直到超过终点。在这种情况下，终点本身就是您应该使用的最小值。
现在，第一步。如果您观察这个方程的一般形状，
y' = -m1·g' / g² = 0
很明显，这只有在
g' = 0的情况下才可能实现。因此，我们必须解决。
g' = 0 =>
m5·a·Xm5-1+ m6·b·Xm6-1+ m7·c·Xm7-1= 0
现在我们在数学上遇到了困难。这个方程有一个微不足道的解 (X=0(t=-m8))。然而，只有在特定参数m5-7的值下才存在一般解决方案。对于那些参数的所有其他值，不存在一般解决方案，您必须通过数值方法找到解决方案。
因此，我建议您使用：
对于第三步，您需要再次重复取导的过程；
y' = -m1·g' / g²
=> y'' = m1·( 2·(g')² - g²·g'' ) / (g²)²
在这里
g'' = (m5-1)·m5·a·Xm5-2+ (m6-1)·m6·b·Xm6-2+ (m7-1)·m7·c·Xm7-2。注意，您只需要知道它的符号，因此只需使用分子即可。
sgn(y'') = sgn( m1·( 2·(g')² - g²·g'' ) )
这是如何在MATLAB中实现的：
g = @(t) 1 + [a b c]*X(t).^[m5; m6; m7];
gp = @(t) ([m5 m6 m7].*[a b c]) * X(t).^[m5-1; m6-1; m7-1];
gpp = @(t) (([m5 m6 m7]-1).*[m5 m6 m7].*[a b c]) * X(t).^[m5-1; m6-1; m7-1];
sgn = @(t) m1*(2*(gp(t)).^2 - g(t).^2*gpp(t) );
然后，使用之前找到的所有根
R，您只需评估所有内容：
extrema = y(R);
signs = sgn(R);
y0 = y(t0);
现在，您已经获得了所有根，所有在根和端点处的
y值以及二阶导数的所有符号....
目标2已完成！
t在哪里？ - Dan