Weird behavior regarding pi in symbolic expressions
58 views (last 30 days)
yetanothernickname on 27 Dec 2018
I have come across some weird behavior where I'm curious whether somebody might know the cause.
I'm trying to check if two (rather complicated) symbolic expressions are the same (spoiler: they are) using a matlab live script
They both contain pi as a part of a fraction. Strangely in one of them pi internally gets replaced by a numeric representation, as can be seen in the output. (More specifically 2/pi^2 gets replaced by 562949953421312 / 2778046668940015)
Thus when I substract the two expressions from each other and use simplify() on the result its not zero (as it should be).
This can be fixed, by defining pi as my own symbol, by adding the line "syms pi". This prevents Matlab from replacing pi and the result of the final substraction is indeed 0.
Does anybody know the cause of this odd behavior? Why does Matlab keep pi in one expression, but unnecessarily replaces it in the other?
This is the code where the behavior can be observed (at least for me, as mentioned I'm running it as a .mlx live scrpt)
(You can ignore most of the stuff in the front, the thing that matters is that rest is !=0 in the first run, but after adding "syms pi" rest=0)
syms phi phid phidd ra ri b l mb ml J g di Dra Dri Db
syms Mpressure p0 lr rho eta lambda Re deltapi deltapr
syms A r d V0
syms phiddold phiddnew Mpressold Mpressnew
syms Reold Renew pralt prnew
A = b*(ra-ri)
d = (ra - ri)
r = (ri + (ra-ri)/2)
Dra = 0
Dri = 0
Db = 0
V0 = pi*(ra^2-ri^2)*b
phiddold = (-(A*r*(deltapr - p0) - (l^2*mb*phid^2*sin(phi))/4 - (l^2*ml*phid^2*sin(phi))/8 + g*l*mb*cos(phi/2) + g*l*ml*cos(phi/2))/(J/2 + (l^2*mb)/2 + (3*l^2*ml)/8 + (l^2*mb*cos(phi))/2 + (l^2*ml*cos(phi))/4 + (2*A*V0*lr*r*rho)/(di^2*pi^2)))
phiddnew = -(g*l*mb*cos(phi/2) - (l^2*ml*phid^2*sin(phi))/8 - (l^2*mb*phid^2*sin(phi))/4 + g*l*ml*cos(phi/2) + (b*(ra - ri)*(ra/2 + ri/2)*(pi*deltapr*di^2 - pi*di^2*p0 + 2*lr*phid^2*rho*Db*ra^2 - 2*lr*phid^2*rho*Db*ri^2 + 4*lr*phid^2*rho*Dra*b*ra - 4*lr*phid^2*rho*Dri*b*ri))/(di^2*pi))/(J/2 + (l^2*mb)/2 + (3*l^2*ml)/8 + (l^2*mb*cos(phi))/2 + (l^2*ml*cos(phi))/4 + (b*(2*lr*rho*b*ra^2 - 2*lr*rho*b*ri^2)*(ra - ri)*(ra/2 + ri/2))/(di^2*pi))
rest = simplify(phiddold-phiddnew)
Walter Roberson on 27 Dec 2018
When you do not use
then when you use pi in an expression, it is the numeric approximation to pi that is used, 3.141592653589793115997963468544185161590576171875 .
When you construct symbolic expressions, the parts that are numeric are constructed as numeric in context. Then at the point where the numeric part would be combined with a symbolic value, the numeric point is converted to symbolic. The process of conversion to symbolic examines some special values such as 3.141592653589793115997963468544185161590576171875 and simple ratios of that, and recognizes them specially and converts them to special expressions. The ones that do not match the special patterns of expressions, are put through a conversion that attempts to match them as expressions in square roots of rationals and rationals, to within tolerance.
In the first of the two expressions, the (di^2*pi^2) is converted to (numeric pi)^2 * di^2 and then the (numeric pi)^2 is not recognized as special and so is converted to the rational 2778046668940015/281474976710656
In the second expression you have (di^2*pi) which is converted to (numeric pi) * di^2 and then the (numeric pi) is recognized as special and converted to symbolic pi.
So... assume that pi will be treated numerically in symbolic expressions, but that from time to time it might be replaced by the symbol.
And if that is not what you want, then assign the meaning you want to Pi and use that instead. For example,
Pi = sym('pi'); %if you want it treated as symbolic
Pi = sym(pi, 'f') %if you want it treated as 884279719003555/281474976710656
Pi = sym(pi, 'd') %if you want it treated as symbolic 3.1415926535897931159979634685442
but in general remember that any numeric value in a symbolic expression might get converted to symbolic form, like sym(sqrt(2)) will be converted to sym(2)^(sym(1)/sym(2)) . You should avoid mixing floating point expressions like sqrt(2) with symbolic expressions and should force treatment of numeric constants according to the way you want.