Weird behavior regarding pi in symbolic expressions

Hey all,

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)

clear
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
%syms pi

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)

[Stephen23]

"To create a symbolic number representing the constant π, use a = sym(pi) instead."

https://www.mathworks.com/help/symbolic/sym.html#mw_8232aaac-d9f8-4d4d-a4c1-4779771a2d61

Answer 1

[Walter Roberson]

When you do not use

syms pi

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.

[yetanothernickname]

Thank you for the great explanation!

[Hank]

Yes, thanks

Is there a way to force some expression containing pi to be interpreted as part of the symbolic part of the equation? It seems a shame to have to define my own Pi just because i'm using some term like (pi/2)^(-1/4).

Is there some clever way to write:

syms E(x) w
E(x) = (pi/2*w^2)^(-1/4) * exp( -(x/w)^2 )

so that the entire expression is interpreted symbolicly?

[Hank]

To answer my own question, I found a fix on another forum.

Unfortuneately its still a a little ugly but we can say:

syms E(x) w
E(x) = (sym(1)*pi/2*w^2)^(-1/4) * exp( -(x/w)^2 )

The sym(1) term creates a unity symbolic object in the expression and forces the whole thing to remain symbolic.

Cheers, H

[Walter Roberson]

I typically

syms pi

If I am going to be using pi symbolically.

Answer 2

[madhan ravi]

Usually to convert symbolic to double

vpa() % is a workaround

Answer 3

[Alec Jacobson]

@Walter Roberson's

syms pi
Pi = sym('pi');

answers aren't 100% satisfying because they are not really associating that variable with the symbolic value of π, rather just making an arbitrary variable named pi. This is verified by trying:

syms pi
cos(pi)

which doesn't simplify.

Instead if you use something like

Pi = str2sym('acos(-1)')

you get a variable Pi that refers to the value π symbolically. This can be verified:

Pi = str2sym('acos(-1)')
cos(Pi)

which returns

ans =

-1

(as a symbolic, i.e., not float, value).

The other issue with

syms pi

Is that it clobbers the built-in numeric value for pi. So any time you build a function from your symbolic expression it will be confused.

syms pi
f = matlabFunction(cos(pi*x))

you get a function of x and this "variable" pi:

f =

function_handle with value:

@(pi,x)cos(pi.*x)

which you'll quickly realize you can't even call with `f(pi,x)` because pi still refers to the symbolic variable.

Meanwhile,

Pi = str2sym('acos(-1)')
matlabFunction(cos(Pi*x))

produces the expected numerical function:

ans =

function_handle with value:

@(x)cos(x.*pi)

str2sym('acos(-1)') is really awkward. I'm not sure if there's a better way to get at the symbolic value π. Perhaps someone with more experience with the symbolic toolkit knows.

[Alec Jacobson]

Ah, seems as long as you haven't already made pi a symbolic variable (e.g., using syms pi) then you can use pi as a string. Probably a safe way is to do:

clear pi
Pi = str2sym('pi');

[Walter Roberson]

Pi = sym(pi);
Pie = sym('pi')

Pie = 

π

syms pi theta
pi+theta

ans = 

vpa(ans)

ans = 

Pi+theta

ans = 

vpa(ans)

ans = 

Pie + theta

ans = 

vpa(ans)

ans = 

Pie - pi

ans = 

0

So when you syms pi or you assign a value sym('pi') or str2sym('pi') then you get a non-special variable that happens to display as π here but is not the transcendental constant.

If you sym(pi) while pi is still the floating point value, or you sym() a floating point value within roughly 6e-14 of that value, then the symbolic engine converts it to a reference to the transcendental constant that the symbolic engine knows how to reason about.

[Alec Jacobson]

interesting that `sym(pi)` does that conversion. I guess it's sadly not applied on the symbol level in an expression involving the built-in pi and a symbolic variable (like in the original post or) here:

clear pi
x = sym('x')
sqrt(pi)*x

produces an inexact result

ans =

(3991211251234741*x)/2251799813685248

but creating a variable for the constant

Pi = str2sym('acos(-1)');
x = sym('x');
sqrt(Pi)*x

produces the exact result

ans =

x*pi^(1/2)

[Walter Roberson]

sqrt(pi)*x has the sqrt(pi) being evaluated at double precision, returning a double-precision value. The symbolic engine is not asked to convert that double precision number until it is asked to execute times(NUMERIC_SQRT_PI, x) . And by that time, NUMERIC_SQRT_PI is too far away from the patterns that the symbolic engine knows to look for.

Every double precision value could potentially be expressed as a algebraic number in relationship to pi. 2.01659265358979 could potentially be reconstructed as pi - 9/8 . sqrt(2) could potentially be expressed as 562949953421312/1250560371546297 * pi .

Obviously, it would not be a good idea to convert every number in relationship to pi . It might not be unreasonable to expand the search to include "simple" rational numbers times sqrt(pi) or pi^2, but it is not immediately clear to me where the hard boundary should be... why not include pi^3 or 4th root of pi or ... ? Is there a branch of statistics or physics that commonly uses sqrt(pi) in calculations but not 4th root of pi ?

[Stephen23]

x = sym('x')

x = 

x

p = sym(pi)

p = 

π

y = sqrt(p)*x

y = 

z = subs(y,x,1)

z = 

double(z)

ans = 1.7725

Answer 4

[Steven Lord]

See the description of the flag input argument on the documentation page for the sym function. The default method to convert a number into a sym object (the 'r' flag) recognizes certain forms of numbers, of which a rational multiple of pi is one, and represents that not necessarily as the exact floating point number but as the form that was used to compute it. From that documentation: "When sym uses the rational mode, it converts floating-point numbers obtained by evaluating expressions of the form p/q, p*pi/q, sqrt(p), 2^q, and 10^q (for modest-sized integers p and q) to the corresponding symbolic form."

So if you convert something like 3*pi/4 to a symbolic value, it recognizes that second form (with p = 3 and q = 4, both "modest-sized integers") and returns .

s = sym(3*pi/4)

s = 

Something like pi^2 is not one of the forms that gets recognized so instead the 'f' conversion method is used, which "returns the symbolic form for all values in the form N*2^e or -N*2^e, where N >= 0 is a nonnegative integer and e is an integer. The returned symbolic number is a precise rational number that is equal to the floating-point value. For example, sym(1/10,'f') returns 3602879701896397/36028797018963968."

s = sym(pi^2)

s = 

[numS, denS] = numden(s)

numS = 

2778046668940015

denS = 

281474976710656

Using the terminology in the description of the 'f' flag, e is -48.

e = -log2(denS)

e = 

If we compute a rational approximation to pi^2 in double precision, we should find that the numerator and denominator of s are the same multiple of the numerator and denominator of that rational approximation, respectively.

[num, den] = rat(pi^2, eps)

num = 107607675

den = 10902937

vpa([numS/num; denS/den])

ans = 

I'd say that's close enough.

[Walter Roberson]

By the way,

sym([306/133*pi, 307/133*pi, 308/133*pi, 307/132*pi, 307/134*pi])

ans = 

That is, 307/133*pi gets converted to a rational rather than a multiple of pi. 133 as the denominator is the smallest rational I found that does not convert to a multiple of pi. on the flip side, rationals less than 1, the first I find is

sym(91/509*pi)

ans =