MATLAB Examples

Polysym User Guide

The polysym class is a symbolic manipulation package that helps users to put equations in a form that Bertini can read. It is similar in behavior to the sym class in the Symbolic Math Toolbox, but limited to commands that are compatible with polynomial systems. Since its purpose is to create strings for Bertini to read, it has limited capabilities compared to the Symbolic Math Toolbox; but it can perform several kinds of arithmetic and matrix operations.

Contents

Construction

A polysym object can be created using a command like:

x = polysym('x') %#ok<*NASGU,*NOPTS> <-- ignore the green stuff (editor options)
x = 

x

This creates an object x which is displayed as 'x'. An array of terms with subscripts can be created by adding a second argument for the dimensions:

v = polysym('x',[1 3])
v = 

x1 x2 x3

If a scalar is provided for the dimension, a square matrix is created:

x = polysym('x',2)
x = 

x1_1 x1_2
         
x2_1 x2_2

A shortcut can be used to define multiple scalars, equivalent to x=polysym('x'), etc.:

polysyms x y z

A variety of inputs are accepted by polysym:

x = polysym('x');                   % char
x = polysym(pi);                    % numeric
x = polysym(polysym('x'));          % polysym
x = polysym(sym('x'));              % sym
v = polysym({'x', 0, vpa(pi)});     % cell array of any of the above types

Arithmetic operations

Polysym objects can be combined using arithmetic operators:

v = polysym('v',[1 2]); w = polysym('w',[1 2]);
disp(v+w)
disp(v-w)
disp(v./w) % element bw element operations
disp(v.*w)
v1+w1 v2+w2
v1-w1 v2-w2
v1/w1 v2/w2
v1*w1 v2*w2

(Of course, v./w does not result in a polyomial expression; division will be used mainly to calculate numerical coefficients.) The above operations can also involve numbers, e.g.,

disp(v*2)
disp(v.^2)
v1*2 v2*2
v1^2 v2^2

Matrix multiplication is allowed, but not matrix division.

M = polysym('v',2); L = polysym('w',2);
z = M*L;
z = M^2;
try
    z = M/L;
catch err
    disp(err.message)
end
Right division (/) must have scalar divisor.

Other arithmetic operations implemented for polysym objects are sum, cumsum, prod, cumprod, and diff. Also overloaded are cos, sin and exp:

y = exp(x)
y = 

exp(x)

Conversion methods

The constructor for polysym is able to convert strings, numbers, objects from the Symbolic Toolbox and cell arrays into the polysym class. There are also methods to convert them back: cellstr, double, sym, vpa:

polysyms x y z
disp(cellstr([x y z]))
disp(sym([x y z]))
    'x'    'y'    'z'

[ x, y, z]
 

Implied conversions can occur when a polysym object is assigned to a member of some other class. In the following example, the polysym object is implicitly converted to double:

M = zeros(2);
M(1) = polysym('1');
disp(M)
     1     0
     0     0

Variable-precision numbers can be converted to polysym and back with little loss of precision, as long as the precision is not changed:

digits(32);         % set precision to 32 significant figures
x = vpa(pi);
y = vpa(polysym(x));
disp([x; y])
 3.1415926535897932384626433832795
 3.1415926535897932384626433832795
 

They are not exactly the same, but they agree to the expected precision, which is 32 significant figures:

disp(x-y)
2.8841974850346244501170641905768e-33
 

The number of digits of precision can be obtained or changed using the command digits. Also, as in the command vpa, the user can specify a different precision for a single number using the command vpa(p,d). However, this can lead to unexpected loss of precision as any conversions are done in the precision determined by digits:

x = vpa(pi,40); % use 40 digits for pi
disp(x)
disp([x; x(1)])
3.141592653589793238462643383279502884197
 
 3.1415926535897932384626433832795
 3.1415926535897932384626433832795
 

In the above example, one would expect x and x(1) to be the same, but x(1) is displayed to the default 32 digits instead of 40. The stored precision is actually the same:

disp(x-x(1))
0.0
 

However, real loss of precision can occur if x is converted to polysym and back. Below, the difference between x and y is much greater than the expected one part in $10^{40}$:

p = polysym(x);
y = vpa(p);
disp(y-x)
-2.8841971640237928394121498415647e-33
 

Thus, it is generally better to set the desired precision in advance instead of specifying it for individual numbers.

Other methods

Some other methods that are overloaded and work much like their builtin counterparts:

  1. relational operations (==, ~= and isequal) are allowed, but not > or <.
  2. state identification: isempty, isone, iszero
  3. array creation and concatenation operations: diag, cat, horzcat and vertcat
  4. matrix operations: cross, tril, triu, transpose (.'), trace
  5. sorting and reshaping arrays: blkdiag, sort

These have very similar behavior to their builtin counterparts. Many other builtin operations such as unique work on polysym objects without overloading.

Constants

The square root of negative 1 is represented as I and pi is represented as Pi in Bertini. Equivalent expressions are converted:

disp(polysym(1i))
disp(polysym(3+4i))
disp(polysym(sqrt(-1)))
disp(polysym('pi'))
disp(polysym(pi))
1*I
3+4*I
1*I
Pi
3.14159265358979

Note that in the last example, pi is the MATLAB builtin constant, so its numerical value is passed to polysym.

Formatting

Finally, polysym has methods for formatting expressions according to the rules of syntax in Bertini input: format_array to format an array for display, with options for the delimiter between terms and the justification; and [format_equations| to format equations.

p = polysym('x',2);
disp(p.format_array)
x1_1 x1_2
x2_1 x2_2
f = polysym('f',[4 1]);
disp(format_equations(f,p(:)))
    'f1=x1_1'
    'f2=x2_1'
    'f3=x1_2'
    'f4=x2_2'

disp(p.format_array(' ','left'))
x1_1 x1_2
x2_1 x2_2

Use in functions

One of the main purposes for polysym is to allow the user to create equations for input into Bertini using functions that also work for ordinary numbers. For example, the following function gives the expansion of exp(x) in x:

f = @(x,n) sum(x.^(0:n)./factorial(0:n));
numeric_value = f(1,10)
Taylor_expansion = f(polysym('x'),4)
numeric_value =

    2.7183


Taylor_expansion = 

1+x+x^2/2+x^3/6+x^4/24

Pitfalls to avoid

A little care must be used in the function file to avoid pitfalls:

1. As is true for sym objects, polysym objects cannot be used with relational operators such as the greater than (>) sign.

try
    polysym(5) > polysym(4); %#ok<VUNUS>
catch err
    disp(err.message)
end
Undefined function 'gt' for input arguments of type 'polysym'.

2. Arrays that will be assigned polysym objects cannot be pre-allocated using commands like zeros because changes to this array would involve attempts to assign a POLYSYM object to a numerical array.

v = zeros(2);
disp('x(1) = polysym(''1'')')
try
    v(1) = polysym('1');
catch err
    disp(err.message)
end
x(1) = polysym('1')
If possible, multiply a |polysym| object by zero instead. For example,
the function below can take inputs from a variety of classes and output
an array of zeros of the correct class:
makeZeroArray = @(x) 0*x;
y = makeZeroArray(1:5)
y = makeZeroArray(polysym(1:5))
disp(['class(y) = ',class(y)])
y =

     0     0     0     0     0


y = 

0 0 0 0 0
class(y) = polysym

3. MATLAB evaluates expressions from left to right. In an example below, a polysym object has a fractional coefficient. If the fraction is to the left of the polysym object, it is evaluated to double precision. To get the exact symbolic representation of the fraction, put the polysym object on the left.

v = polysym('x');
disp('9/16*x = ')
disp(9/16*v)
disp('x*9/16 = ')
disp(v*9/16)
9/16*x = 
0.5625*x
x*9/16 = 
x*9/16

4. The polysym class does not support complex conjugation, since it does not give rise to polynomial expressions in complex variables.

Since many users, accustomed to dealing mainly with real numbers, use the symbols for a complex conjugate transpose (') and non-conjugate transpose (.') interchangeably, polysym will accept the operator (') but return the non-conjugate transpose along with a warning:

p = polysym('x',2);
q = p';
Warning: You tried to calculate a complex conjugate transpose (e.g., p'). The
POLYSYM class does not support complex conjugation, so the non-conjugate
transpose (.') was substituted. If this is correct, consider changing (') to
(.') in your code. 
isequal(q,p.')
ans =

     1

Functions such as dot and abs that use complex conjugation are not supported.