Weird behavior regarding pi in symbolic expressions

115 次查看（过去 30 天）
yetanothernickname
yetanothernickname 2018-12-27
评论： Walter Roberson 2023-7-19
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)

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2018-12-27
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.
  4 个评论
显示 2更早的评论
Hank
Hank 2019-7-11
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

请先登录，再进行评论。

更多回答（3 个）

madhan ravi
madhan ravi 2018-12-27
Usually to convert symbolic to double
vpa() % is a workaround
Alec Jacobson
Alec Jacobson 2023-7-18
@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.
  4 个评论
显示 2更早的评论
Alec Jacobson
Alec Jacobson 2023-7-19
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
Walter Roberson 2023-7-19
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 ?

请先登录，再进行评论。

Steven Lord
Steven Lord 2023-7-19
Ran in:
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.
  1 个评论
Walter Roberson
Walter Roberson 2023-7-19
Ran in:
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 = 

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Variables, Expressions, Functions, and Preferences 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by