# 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 :

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:

- relational operations (
`==`,`~=`and`isequal`) are allowed, but not`>`or`<`. - state identification:
`isempty`,`isone`,`iszero` - array creation and concatenation operations:
`diag`,`cat`,`horzcat`and`vertcat` - matrix operations:
`cross`,`tril`,`triu`,`transpose (.')`,`trace` - 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.