matlab for automatic differentiation

Question

if someone of you could help me out or point me in the right direction that would be wonderfull.

I have the following formula which has 8 parameters that differ for each subject. The equation describes a growth pattern with t being time and y being the growth. m1 to m8 are the parameters that differ for each subject

y = m1*(1-1/(1+(m2*(t+m8))^m5+(m3*(t+m8))^m6+(m4*(t+m8))^m7))

What i would like to do is insert this formula and calculate (numerically) the velocity curve. From this information i would like to obtain the following

• the maximum velocity and the corresponding age and height at that time
• the minimum velocity before the maximum velocity is reached.

and this in such a way that i just give a list of all my subjects and their parameters and just receive my results?

I have no background at all in programming and mathematics but i would like to be able to do this. I have matlab to my disposal. I'm trying to figure out how things work but i can't grip my head around it. As a biologist I will follow a programming course next year ;-)

Can anyone of you help me out.

Thanks

Update but still not there (thanks to you guys :-))

-I have made a tab delimited file derived from excel and having the following structure column A: the ID of the subject Colum B to I: the values of the different parameters m1 ... m8 (as in the equation) each row is a different subject with different parameters

• I have chosen import selection and next generate function (or should i choose generate script).

• If i open this i see a new numbered tab appear. First line is afunction[ID,m1,...,m8]=importfile. then follow 112 lines.

• I have copied the text (thanks to dan) as given by dan starting from line 113.

  tmin=0; tmax=20; dt=1/12; t=timn:dt:tmax; y = m1.(1-1./(1+(m2.(t+m8)).^m5+(m3.(t+m8)).^m6+(m4.(t+m8)).^m7)); dy=diff(y)./dt; max(dy); min(dy); imax=find(dy==max(dy))+1; imin=find(dy==min(dy))+1; t(imax); t(imin); y(imax); y(imin);

  • next i clicked run and was asked to save the file wich i've done.
  • now i clicked run again but i get an error

Anyone that can point me in the right direction?? Thanks a lot

Greetings

Answer 1

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:

1. Find the velocity curve
2. Find the maximum y(tmax) = ymax
3. Find the minimum ymin, tmin that occurs before the maximum ymax

When translated into mathematical lingo:

1. Find the derivative y' = dy/dt
2. (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:

1. find all the roots of y' = 0.
2. evaluate y on these roots
3. evalate y'' on these roots
4. 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

• roots when m5-7 are integers
• FindRealRoots when they are not

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!

Answer 2

well as your question is a bit general I can just provide some pointers where to look.

• for the formula you could store in a [function][1] or [function handle][2] - also there is a symbolic math toolbox that you could use (in that case you could use symbolic differentiation to find the extrema)

• Next is basically loops (see for or while) - will depend a bit in how you read your data into matlab.

• last would be finding the extrema (I suppose you don't know the exact time of them) - lookup the fmin command for numeric minimization (is for minima - but for maxima just use -your function as objective) - if there are several you might need to play a bit with it's settings.

• lets not forget some output. I think xlswrite is most comfortable. But help will point you to alternatives.

---

as a side note you should try to do the differentiation manually to check how many solutions you should expect.

Answer 3

If you equation is of the form: y = f(t,...,... ...) and all the other parameters are held constant over all time for each subject then you can do it like this.

First decide on a time range (so I'm going with 10 years, but that's totally dependent on your work):

tmin = 0;
tmax = 10;

Then you must choose a time resolution (I've gone with monthly, but you'll probably want something finer):

dt = 1/12;

Now create a t vector and solve for y numerically:

t = timn:dt:tmax;
y = m1.*(1-1./(1+(m2.*(t+m8)).^m5+(m3.*(t+m8)).^m6+(m4.*(t+m8)).^m7));

First order numerical differentiation is just finding the slope between each two consecutive points of y. I.E. (y(i) - y(i - 1)) / (t(i) - t(i - 1)) and note that t(i) - t(i - 1) is always equal to dt. So in matlab we can just do this:

dy = diff(y)./dt;

but now dy has one less element than dt does. Keep this in mind.

So to find the max and min velocity:

max(dy);
min(dy);

and to find the corresponding times and heights:

imax = find(dy==max(dy)) + 1;
imin = find(dy==min(dy)) + 1;

t(imax);
t(imin);
y(imax);
y(imin);