Thread Subject:

Symbolic expression problem

Help is needed! I have a problem with symbolic derivative expression by Matlab. I work with Hill equation. When I evaluate the symbolic expression of its first derivative, I found out the value isn't quite correct. I diff() with respect to the correct variable, vectorize the expression (If calculated correctly, c=10 should correspond to half amplitude, also c=10 is where the derivative is at max). I don't know what else is there to be done. The following code is what I have done:

%Hill's equation
%
% c Independent variable.
% q Hill's coefficient.
% s^q Half-saturation term.
% Y Amplitude

% Y = (c^q)/(c^q + s^q) ---This is the pseudocode for Hill equation----


clear all;

c=(0:1:100)';
q=2;
s=10;

y = (c.^(q))./((c.^q)+s^q);

dy_dx = y./c; %Calculate slope at every point. This is gives correct result.

dy_dx_symtb = c.^q.*q./c./(c.^q+s.^q)-(c.^q).^2./(c.^q+s.^q).^2.*q./c;
%This is generated by Matlab symbolic toolbox
% by the following commands: ans=diff(y, c); vectorize(ans).

%Plot results
figure(1); subplot(211); plot(c, y); xlabel('c'); ylabel('y'); legend('Hill coefficient = 2');
figure(1); subplot(212); plot(c, dy_dx, c, dy_dx_symtb); xlabel('c'); ylabel('y'); legend('calculated dy/dx is correct', 'evaluating symbolically derived dy/dx');

KC wrote:
> %Hill's equation
> %
> % c Independent variable.
> % q Hill's coefficient.
> % s^q Half-saturation term.
> % Y Amplitude

> % Y = (c^q)/(c^q + s^q) ---This is the pseudocode for Hill equation----

> c=(0:1:100)';
> q=2;
> s=10;

> y = (c.^(q))./((c.^q)+s^q);

> dy_dx = y./c; %Calculate slope at every point. This is gives correct result.

That step is obviously incorrect.

Consider the first c value, c = 0; y at c=0 is going to be 0^q on the numerator and
0^q+s^q on the denominator, and as long as s is non-zero and q is non-zero, that's going
to be 0 on the numerator and non-zero on the denominator, which is going to give
a well-defined 0 value. With q=2, the numerator is everywhere smooth, and with non-zero
s on the denominator, the denominator can generate a discontinuity only if q is odd
or c or s are imaginary numbers, neither of which condition holds in this instance.
Thus we see that the derivative must be smooth for Y when q=2 and s and c real numbers.

But consider your dy_dx at c=0: you have 0 on the numerator and 0 on the denominator,
which is going to give NaN in matlab (and would give an error if it were evaluated
at the symbolic level.) So we can see by inspection that your dy_dx is incorrect
at at least one point, which leads us to suspect the results.

You cannot calculate dy_dx as y / x: you have to calculate it as
(y(t+delta) - y(t)) / (t + delta - t)
Your delta would be 1 in the case of 0:1:100 so your dy_dx should
be (diff(y) ./ 1) rather than y / c.