Obtaining accurate derivatives using Symbolic Math Toolbox

Hi all

I'm coding a model that requires extensive partial differentiation (first and second-order derivatives). I decided to try using Symbolic Math Toolbox for the partial derivatives. I used the following methodology:

1) Use the diff() function to determine derivatives

2) Use subs() to determine the value of the derivative at a certain point

3) Use double() to convert symbolic values to numerical ones that can be used in subsequent calculations

In order to verify that round-off error or internal cancellation don't reduce the accuracy of my derivatives obtained with steps 1 -3, I checked my results against numerical derivatives (formulae for numerical first and second partial derivatives w.r.t. x and y provided at https://en.wikipedia.org/wiki/Finite_difference under the heading "Multivariate finite differences". h=k=0.00001 used).

I found that the % difference between some of the "analytical" derivatives obtained using Symbolic Math Toolbox (steps 1-3) and the checks (numerical derivatives) is at times not satisfactory (up to 0.0076%, would be preferable if this difference were at least 2 orders of magnitude smaller).

How do I increase the agreement between the numerical and the "analytical" derivatives?

I have tried using digits(300) to increase the precision used in calculations. I've also used vpa(...,300) to prevent truncation errors during differentiation.

Is it a fool's errand to try to get the numerical and the "analytical" derivatives to match up to a greater extent?

Thanks in advance.

Alex

Answer 1

[Andreas Apostolatos]

Hello,

The evaluations of symbolic expressions that represent analytic derivatives in MATLAB are only governed by round-off errors that might occur during the evaluations. The 'vpa' functionality can provide you with higher accuracy in terms of the round-off errors of the corresponding evaluations. Therefore, only the round-off errors play a role when evaluating a symbolic expression in MATLAB.

On the contrary, when computing a derivative numerically, e.g., using the finite difference method, then the solution accuracy depends primarily on the discretization parameters, in your case 'h' and 'k', and decreasing those should result in principle into a better approximation of the corresponding derivative. Of course, the round-off errors come into the play in this case as well, and which of the two (discretization vs round-off error) is more significant depends solely on the level of both errors.

Having said that, I would expect that you can get a better agreement than 0.0076% between the analytical and the discretized expression of the underlying derivative, if you would decrease the step sizes 'h' and 'k', if the value of your derivatives is not that high such that round-off errors would be significant at the level of 0.0076%. However, there will certainly be a point when decreasing the step sizes, where the round-off error of the computations is more significant than the discretization error, but I would expect this to occur in a much lower error level than 0.0076%, provided that the actual value of the derivative is not very high.

I hope that this helps with the understanding.

Kind Regards,

Andreas

[Alexander Schulze-Hulbe]

Hi Andreas

Thank you very much for your helpful answer. As you suggested, I decreased 'h' and 'k', such that h=k=1e-6. This produced slighly better results, but h=k=1e-7 made the %error increase. This suggests that the round-off error becomes more significant than the error introduced by the discretization for this value of 'h'/'k'.

You made me realise that the problem may well lie with the numerical derivatives - perhaps comparing Symbolic Math toolbox results with them is not the best way to verify that my model is accurate.

I have also realised that both the output of the numerical as well as the symbolic derivatives obtained with steps 1-3 discussed above sometimes have irrational parts. Could you perhaps explain why that is (should I rather ask this in a separate post)?

Kind regards

Alex

[Andreas Apostolatos]

Hello Alex,

I am happy to have helped with the understanding!

Please keep in mind, that decreasing the step size in the frame of the finite difference method will eventually lead into having a solution dominated by the round-off error, see the following link for more information,

https://en.wikipedia.org/wiki/Numerical_differentiation#Step_size

In fact, there is a threshold for the step size below which the approximation of the derivative will deteriorate. Thus, one may speak of an optimal step size in the sense that it provides the best possible approximation of the derivative. There are different finite difference formulations that you can try, each of which coming with its own benefits and shortcomings.

I think that comparing the numerical approximation you obtain by means of the finite difference method against the analytic solution that you obtain against the Symbolic Math Toolbox is productive. The analytic solution can be safely assumed as a reference solution in that sense. The issue is that the finite difference scheme provides a convergent approximation to the reference solution up to a certain level of the step size, as mentioned above.

You mentioned that sometimes you obtain irrational parts in the solution from both the analytic derivation by means of the Symbolic Math Toolbox and the numerical solution. I assume that by irrational parts you mean such numbers as 'sqrt(2)'. I do not immediately see any restriction on having such values in the analytic derivative of a general function. Regarding the numerical expression, it might be that the Symbolic Math Toolbox simplifies the underlying expressions leading to various parts including irrational parts. Therefore, I would not exclude this possibility.

I hope this helps as well!

Kind Regards,

Andreas

[Alexander Schulze-Hulbe]

Hi Andreas

Thank you very much - this was also very helpful!

Apologies - I meant to say imaginary parts in the solution in my previous comment, not irrational ones - why does MATLAB sometimes give me answers with non-zero imaginary components (e.g. 5 +0.1*i) for both numerical and analytic derivatives?

Kind regards

Alex

[Andreas Apostolatos]

Hi Alex,

Thanks for the clarification. It seems that some MATLAB functions may result in complex expressions while being differentiated due to automatic simplifications by means of standard rules, see the following MATLAB Answers thread accordingly,

https://www.mathworks.com/matlabcentral/answers/398706-derivative-of-a-real-function-becomes-complex

There it is also suggested a workaround for a particular case, where the square root of an expression is converted into an expression containing an absolute value by means of standard trigonometric conversion rules. Thus, in case your expression is simplified leading into an expression that contains imaginary parts, you can enforce your expression not to use the alternative simplification by means of command 'rewrite'. For more information accordingly please visit the following documentation page,

https://www.mathworks.com/help/symbolic/rewrite.html

Another possible workaround is to set assumptions on the symbolic variables at hand. If for instance 'x' is declared as a symbolic variable, namely,

syms x

Then, you can declare variable 'x' as real, by means of the following command,

assume(x,'real')

For more information, please visit the following documentation page,

https://www.mathworks.com/help/symbolic/sym.assume.html

I hope this information makes the matters somewhat clearer.

Kind Regards,

Andreas

[Alexander Schulze-Hulbe]

Hi Andreas,

Thank you very much. That is very informative! In general, it appears that symbolic expressions tend to be complex (i.e. isreal returns zero). The model I'm working with doesn't make physical sense with imaginary values, which is why I need to ensure everything remains real.

Unfortunately, "assume" doesn't seem to do the trick. For example, please see below:

syms m1 x epsilon1 n1 real
EoS = [m1 x epsilon1 n1];
syms T d
assume(d,'real')
isreal(d) %returns true
d = EoS(2).*(1-0.12.*exp(-3.*(EoS(3)).*(1./T))); %define d. 
isreal(d) %returns false

In the example above, is there any way I can avoid d from becoming complex?

Kind regards

Alex

[Andreas Apostolatos]

Hi Alex,

Thanks a lot for sharing this code snippet.

What happens here is the following: Variable 'd' is declared as real, but then it is assigned to an expression which is not necessarily real. Therefore you obtain a complex result when executing line,

d = EoS(2).*(1-0.12.*exp(-3.*(EoS(3)).*(1./T))); %define d. 

Function 'assume' makes assumptions on the variable when it is used in the right-hand side of an expression. If you assign the variable to another expression (by placing it to the left hand-side of an equation), then the assumptions are no more applicable.

In order to achieve the desirable result, please declare the variables used within the latter expression as 'real', namely,

syms m1 x epsilon1 n1 real
assume(x,'real') % declare variable 'x' as real since it is used in the following expression
EoS = [m1 x epsilon1 n1];
syms T d
assume(T,'real') % declare variable 'T' as real since it is used in the following expression
assume(d,'real')
isreal(d) % returns true
d = EoS(2).*(1-0.12.*exp(-3.*(EoS(3)).*(1./T))); % define d. 
isreal(d) % returns true this time

Now since variables 'x' (or equivalently 'EoS(2)') and 'T' are declared as 'real', the expression that you assign to variable 'd' is real.

I hope this information helps you further!

Kind regards,

Andreas

[Alexander Schulze-Hulbe]

Hi Andreas,

Thank you very much! You made me realise a mistake: I forgot to define the variable 'T' as real.

I've realised that I get complex outputs because the functions log and sqrt give complex outputs for symbolic variables, even if those symbolic variables are defined as real. For instance:

syms x y z  real
y = log(x)
isreal(y) %returns false
z = sqrt(x)
isreal(z) %returns false

Is there a way to make y and z real?

Thanks in advance for your time.

Kind regards

Alex

[Andreas Apostolatos]

Hi Alex,

Assuming 'x' to be real is not enough for having expressions 'log(x)' and 'sqrt(x)' to be real.

In fact, you need to have variable 'x' be strictly positive in order to ensure that the expressions 'log(x)' and 'sqrt(x)' are real numbers. The following code snippet should achieve what you are seeking for,

>> syms x y z
assume(x, 'positive')
y = log(x)
isreal(y) % now this returns true
z = sqrt(x)
isreal(z) % now also this returns true

I hope this helps.

Kind Regards,

Andreas

[Alexander Schulze-Hulbe]

Hi Andreas

Oh yes, of course! That was helpful, thanks.

What is giving me trouble is an expression of the kind "log(1-x)" - I can't seem to convince MATLAB that "1-x" should be a positive value. If I'm not mistaken, MATLAB accepts that "1-x" is real, just not that it's positive.

Please see my code snippet below.

clear all %clear all symbolic variables
syms m1 d_sigma epsilon1 n1 N_av T d V

d = d_sigma*(1-0.12*exp(-3*(epsilon1)*(1/T)));

assume(m1,'real')
assume(d_sigma,'real')
assume(epsilon1,'real')
assume(n1,'real')
assume(N_av,'real')
assume(T,'real')
assumeAlso(d,'real')
assume(V,'real')
%the values defined in this paragraph all have physical meaning
% and shouldn't be negative.
assumeAlso(m1,'positive')
assumeAlso(d_sigma,'positive')
assumeAlso(epsilon1,'positive')
assumeAlso(n1,'positive')
assumeAlso(N_av,'positive')
assumeAlso(T,'positive')
assumeAlso(d,'positive')
assumeAlso(V,'positive')

A1 = pi*N_av/(6*V)*n1*m1*(d)^3
B1 = log(A1)
test1 = isreal(B1) %returns true, but I'm actually interested in log(1-A1)

B2 = log(1-A1)
test2 = isreal(B2) %returns false - I think this is so because (1-A2) may be negative,
% meaning that complex numbers cannot be ruled out when computing log(1-A1)

assumeAlso(A1<1) %use this to prevent possibilty of taking natural log of a negative number
B3 = log(1-A1)
test3 = isreal(B3) %Still returns false - why?

Thanks in advance for your time.

Kind regards

Alex

P.S. Seeing that the title of this thread may be a little misleading (originally this was about derivatives), I'm going to repost this last question separately as well.

[Andreas Apostolatos]

Hello Alex,

Indeed it is better to have this question posted in a separate thread in order to have each issue addressed in the appropriate channel for the sake of clarity.

Since you are going to repost this question, please provide a link to the new post herein for the sake of completeness.

Thanks.

Kind Regards,

Andreas

[Alexander Schulze-Hulbe]

Hi Andreas

Here's the link.

https://www.mathworks.com/matlabcentral/answers/796617-symbolic-math-toolbox-obtaining-real-values-from-an-expression-of-the-form-log-1-x?s_tid=srchtitle

Kind regards Alex