First, here's a copy-pastable solution. The rest of my answer is the explanation :)
(You'll need FindRealRoots from the FEX to run this).
trange = [tstart tend] % <your range in t>
% define all intermediate functions
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;
% Find the roots of the derivative
R = FindRealRoots(gp, tstart , tend, 2*ceil(max([m5 m6 m7])));
% values at end-points
y0 = y(tstart);
yend = y(tend);
% if we have some roots, separate minima/maxima and sort
if ~isempty(R)
extrema = y(R);
signs = sgn(R);
% Find maximum
[maximum, Imax] = max(extrema(signs < 0));
tmax = R(Imax);
if isempty(tmax) % no maximum in interval; use endpoint
[~,ind] = max([y0, yend]);
tmax = trange(ind);
end
ymax = y(tmax);
% Find minimum *before* maximum
[minima, Imin] = extrema(signs > 0);
tmin = R(Imin);
if isempty(tmin) % no minimum on interval; use endpoint
[~,ind] = min([y0, yend]);
tmin = trange(ind);
end
if any(tmin < tmax)
% make sure to take only one minimum; the one
% just before the maximum
tmin = tmin(find(tmin < tmax, 1, 'last'));
else
[~,ind] = min([y0, yend]); % no minima before the maximum; use end-point
tmin = trange(ind);
end
ymin = y(tmin);
else % no roots found; just order the end-points
[ymax,ind] = max([y0, yend]);
tmax = trange(ind);
[ymin,ind] = min([y0, yend]);
tmin = trange(ind);
end
Now, the explanation.
While the answer by Dan is a nice first step, I feel I must correct it in a few places.
For starters, numerical derivatives are computed as Dan describes them only when
- you have gridded data (the value of
y(t) at a number of distinct points t)
- It is known that the underlying model function
y(t) does not allow for interpolation
And this is simply not true in your case. When you have an explicit function that you can evaluate at any arbitrary point, numerical derivatives are best computed using finite differences, and normally only when explicit derivatives cannot be found analytically.
As an example, the numerical derivative of z = cos(x) at x = π/4 (using central differences with h = 0.001):
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)
Having said that, we move on to the more important point; the mathematics can be progressed way further than Dan did, before you have to resort to numerical methods. Doing so will perhaps be more tedious, but it will improve your understanding of the problem, of MATLAB, mathematics, and it will improve the accuracy and robustness of the outcomes.
For your function y = f(t), you state you have 3 objectives:
- Find the velocity curve
- Find the maximum
y(tmax) = ymax
- Find the minimum
ymin, tmin that occurs before the maximum ymax
When translated into mathematical lingo:
- Find the derivative
y' = dy/dt
- (and 3.) Find all roots
R of y' and evaluate y and y'' at all R, and also y at the initial time t0.
Let's start with objective 1. You have an explicit function of time y = f(t), which is basically a variant form of a rational function,
y(t) = a + f(t)/g(t)
For brevity, let's just drop the (t)-part. That's much easier to write, and it's understood that y, f and g all depend on t. So,
y = a + f/g
You have the advantage that f(t) = <constant> = a = m1, which facilitates things later on. The expression for the derivative of this kind of function takes this form (which is the quotient rule):
y' = (f'·g - f·g') / g²
and, since f(t) = m1 = <constant>, its derivative with respect to time is zero:
y' = -m1·g' / g²
Now, we need to find the derivative of g. The function g is the polynomial
g = 1 + (m2·(t+m8))m5+ (m3·(t+m8))m6+ (m4·(t+m8))m7
Let's define X = t + m8, a = m2m5, b = m3m6 and c = m4m7. With this, g will look a lot less scary:
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
and therefore,
y' = -m1·g' / g² = -m1 · (m5·a·Xm5-1+ m6·b·Xm6-1+ m7·c·Xm7-1) / (1 + a·Xm5+ b·Xm6+ c·Xm7)²
Objective 1 accomplished! You can now evaluate the velocity curve anywhere with full accuracy. Here's how to implement this function in 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;
Now objective 2:
- find all the roots of
y' = 0.
- evaluate
y on these roots
- evalate
y'' on these roots
- evaluate
y on the end-points of the t-interval
Step 3. is necessary to determine whether a root R of y' = 0 describes a maximum or a minimum of y. In case y''(R) < 0, you are dealing with a maximum, in case y'' > 0, a minimum, and y''(R) = 0 you are dealing with an inflection point, which you may safely ignore for this problem.
Step 4. is necessary, because there may be cases where before the maximum, there is no minimum in the function, until way beyond the end point. in that case, the end point itself is the minimum you should use.
Now, step 1. If you look at the general shape of this equation,
y' = -m1·g' / g² = 0
it should be obvious that this is only possible if g' = 0. Therefore, we have to solve
g' = 0 =>
m5·a·Xm5-1+ m6·b·Xm6-1+ m7·c·Xm7-1= 0
And it is only now that we are getting stuck with our mathematical progress. There is one trivial solution to this equation (X = 0 (t = -m8)). However, there is only a general solution for specific values of the parameters m5-7. For all the other values of those parameters, no general solution exists and you'll have to find solutions numerically.
Therefore, I will recommend you to use
For step 3, you repeat the process of taking derivatives one more time;
y' = -m1·g' / g²
=> y'' = m1·( 2·(g')² - g²·g'' ) / (g²)²
where
g'' = (m5-1)·m5·a·Xm5-2+ (m6-1)·m6·b·Xm6-2+ (m7-1)·m7·c·Xm7-2. Note that you just need the sign of this, so only the numerator will suffice),
sgn(y'') = sgn( m1·( 2·(g')² - g²·g'' ) )
Here's how to do that in 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) );
Then, with all roots R found previously, you just evaluate everything:
extrema = y(R);
signs = sgn(R);
y0 = y(t0);
Now that you have all roots, all values of y at the roots and end-points, and all signs of the second derivative as well....
Objective 2 accomplished!