MATLAB Examples

## Contents

```% Sympoly demos ```

## Various ways to create a sympoly

```% A scalar (zero) sympoly z = sympoly; % Scalar sympolys 'x', 'y', 'u', 'v' created in the current workspace sympoly x y u v % A sympoly (identity matrix) array. The numeric element format is % specified by the command window format style. format short g ayuh = sympoly(eye(3)); % Use deal to replicate a sympoly into several [a,b] = deal(sympoly); % Deal can also create a sympoly array S(1:2) = deal(sympoly('x')); % As can repmat R = repmat(sympoly('x'),2,3); whos ```
``` Name Size Bytes Class Attributes A 1x1 684 sympoly B 1x1 448 sympoly F 1x1 2198 sympoly I 1x1 2302 sympoly I1 1x1 2044 sympoly I2 1x1 2302 sympoly Ic 1x1 530 sympoly M 4x4 128 double N 1x1 2302 sympoly Nc 1x1 530 sympoly P_of_x 1x1 1646 sympoly Q 1x1 15424 sympoly R 2x3 1740 sympoly S 1x2 708 sympoly a 1x1 448 sympoly ans 1x1 880 sympoly ayuh 3x3 2496 sympoly b 1x1 448 sympoly f0 1x1 452 sympoly f1 1x1 452 sympoly f2 1x1 452 sympoly f3 1x1 452 sympoly f4 1x1 452 sympoly fA 1x1 452 sympoly fAB 1x1 454 sympoly fAC 1x1 454 sympoly fAD 1x1 454 sympoly fB 1x1 452 sympoly fBC 1x1 454 sympoly fBD 1x1 454 sympoly fC 1x1 452 sympoly fCD 1x1 454 sympoly fD 1x1 452 sympoly fcenter 1x1 462 sympoly fx 1x1 452 sympoly fxx 1x1 454 sympoly fxxx 1x1 456 sympoly fxxy 1x1 456 sympoly fxy 1x1 454 sympoly fxyy 1x1 456 sympoly fy 1x1 452 sympoly fyy 1x1 454 sympoly fyyy 1x1 456 sympoly lambda 1x1 460 sympoly mux 1x1 454 sympoly muy 1x1 454 sympoly p 1x1 734 sympoly p2 1x1 566 sympoly p3 1x1 542 sympoly p4 1x1 566 sympoly p5 1x1 638 sympoly polymean 1x1 748 sympoly polyvar 1x1 1076 sympoly quotient 1x1 542 sympoly remainder 1x1 448 sympoly sx 1x1 452 sympoly sy 1x1 452 sympoly u 1x1 450 sympoly v 1x1 450 sympoly x 1x1 450 sympoly xy 7x2 112 double xyz 11x3 264 double y 1x1 450 sympoly z 1x1 448 sympoly ```

## Arithmetic between sympolys, add, subtract, multiply, divide.

```% add 1 to x 1 + x ```
```ans = 1 + x ```
```% double times a sympoly 2*y ```
```ans = 2*y ```
```% subtraction, and a simple power (x - y)^2 ```
```ans = y^2 - 2*x*y + x^2 ```
```% More complex expressions (x - 2*y)^3/x + sqrt(y^3) ```
```ans = -8*x^-1*y^3 + y^1.5 + 12*y^2 - 6*x*y + x^2 ```

## Synthetic division

```[quotient,remainder] = syndivide(x^2+2*x-1,x+1) ```
```quotient = 1 + x remainder = -2 ```

## Arrays of sympolys

```[x , y ; 1 , x+y] ```
```ans = Sympoly array has size = [2 2] Sympoly array element [1 1] x Sympoly array element [2 1] 1 Sympoly array element [1 2] y Sympoly array element [2 2] y + x ```
```% Arrays of sympolys v = [1 x y x+y] ```
```v = Sympoly array has size = [1 4] Sympoly array element [1 1] 1 Sympoly array element [1 2] x Sympoly array element [1 3] y Sympoly array element [1 4] y + x ```

## matrix multiplication

```A = v*v' B = v'*v ```
```A = 1 + 2*y^2 + 2*x*y + 2*x^2 B = Sympoly array has size = [4 4] Sympoly array element [1 1] 1 Sympoly array element [2 1] x Sympoly array element [3 1] y Sympoly array element [4 1] y + x Sympoly array element [1 2] x Sympoly array element [2 2] x^2 Sympoly array element [3 2] x*y Sympoly array element [4 2] x*y + x^2 Sympoly array element [1 3] y Sympoly array element [2 3] x*y Sympoly array element [3 3] y^2 Sympoly array element [4 3] y^2 + x*y Sympoly array element [1 4] y + x Sympoly array element [2 4] x*y + x^2 Sympoly array element [3 4] y^2 + x*y Sympoly array element [4 4] y^2 + 2*x*y + x^2 ```

## Selective extraction of terms

The second term

```terms(A,2) ```
```ans = 2*y^2 ```
```terms(A,x^2,'extract') ```
```ans = 2*x^2 ```

Delete a term

```p = (1 + x^2 + x^7)^3 terms(p,x^2,'delete') ```
```p = 1 + 3*x^2 + 3*x^4 + x^6 + 3*x^7 + 6*x^9 + 3*x^11 + 3*x^14 + 3*x^16 + x^21 ans = 1 + 3*x^4 + x^6 + 3*x^7 + 6*x^9 + 3*x^11 + 3*x^14 + 3*x^16 + x^21 ```

## Selective deletion of terms

```B = terms(A,x,'extract') ```
```B = 0 ```

Operations on arrays

```sympoly lambda (rand(3) - lambda*eye(3)) ```
```ans = Sympoly array has size = [3 3] Sympoly array element [1 1] 0.96489 - lambda Sympoly array element [2 1] 0.15761 Sympoly array element [3 1] 0.97059 Sympoly array element [1 2] 0.95717 Sympoly array element [2 2] 0.48538 - lambda Sympoly array element [3 2] 0.80028 Sympoly array element [1 3] 0.14189 Sympoly array element [2 3] 0.42176 Sympoly array element [3 3] 0.91574 - lambda ```

## Even eigenvalues, using det, then roots

```roots(det(hilb(4) - lambda*eye(4))) ```
```ans = 1.5002 0.16914 0.0067383 9.6702e-05 ```

## Sum on any dimension

```sum(v,2) ```
```ans = 1 + 2*y + 2*x ```

## And prod

```prod(A(:)) ```
```ans = 1 + 2*y^2 + 2*x*y + 2*x^2 ```

## Orthogonal polynomials from a variety of familes

```% 3rd and 4th order Legendre polynomials p3 = orthpoly(3,'legendre') p4 = orthpoly(4,'legendre') ```
```p3 = -1.5*x + 2.5*x^3 p4 = 0.375 - 3.75*x^2 + 4.375*x^4 ```
```% Orthogonal polynomials are orthogonal over the proper domain defint(p3*p4,'x',[-1,1]) ```
```ans = 0 ```
```% 2nd and 5th order Jacobi polynomials p2 = orthpoly(2,'jacobi',2,3) p5 = orthpoly(5,'jacobi',2,3) ```
```p2 = -1 - 2*x + 9*x^2 p5 = -0.65625 + 7.2188*x + 14.4375*x^2 - 62.5625*x^3 - 31.2812*x^4 + 93.8438*x^5 ```
```% Orthogonal polynomials are orthogonal over the proper domain. % Numerical issures left this just eps shy from zero. defint(p2*p5*(1-x)^2*(1+x)^5,'x',[-1,1]) ```
```ans = -1.4211e-14 ```

## Roots of the derivative of a sympoly

```sort(roots(diff(orthpoly(6,'cheby2')))) ```
```ans = -0.79821 -0.44293 0 0.44293 0.79821 ```

## Error propagation through a sympoly

```% Given a unit Normal N(0,1) random variable, compute the % mean and variance of p(x) = 3*x + 2*x^2 - x^3 sympoly x [polymean, polyvar] = polyerrorprop(3*x + 2*x^2 - x^3,'x',0,1) ```
```polymean = 2 polyvar = 14 ```
```% Compute the mean and variance of x*y + 3*y^3, where x and y are % respectively N(mux,sx^2), and N(muy,sy^2) sympoly x y mux muy sx sy [polymean,polyvar] = polyerrorprop(x*y+3*y^3,{'x' 'y'},[mux,muy],[sx,sy]) ```
```polymean = 9*muy*sy^2 + 3*muy^3 + mux*muy polyvar = 135*sy^6 + sx^2*sy^2 + 324*muy^2*sy^4 + muy^2*sx^2 + 81*muy^4*sy^2 + 18*mux*sy^4 + 18*mux*muy^2*sy^2 + mux^2*sy^2 ```

## A simple construction for a Newton-Cotes integration rule

Here I'll generate Simpson's 3/8 rule.

Simpson's 3/8 rule

```M = vander(0:3); sympoly x f0 f1 f2 f3 % an interpolating polynomial on this set of points % { (0,f0), (1,f1), (2,f2), (3,f3) } P_of_x = [x^3, x^2, x, 1]*pinv(M)*[f0;f1;f2;f3]; ```

## sympoly uses the command window format style to write out the coefficients

Here, I'll force it to be in rational form

```format rat ```

integrate the polynomial over its support

```defint(P_of_x,'x',[0 3]) ```
```ans = 3/8*f3 + 9/8*f2 + 9/8*f1 + 3/8*f0 ```

Or here, a 4 point open Newton-Cotes rule

```M = vander(1:4); sympoly x f1 f2 f3 f4 % an interpolating polynomial on this set of points % { (1,f1), (2,f2), (3,f3) (4,f4) } P_of_x = [x^3, x^2, x, 1]*pinv(M)*[f1;f2;f3;f4] ```
```P_of_x = -1*f4 + 11/6*f4*x - 1*f4*x^2 + 1/6*f4*x^3 + 4*f3 - 7*f3*x + 7/2*f3*x^2 - 1/2*f3*x^3 - 6*f2 + 19/2*f2*x - 4*f2*x^2 + 1/2*f2*x^3 + 4*f1 - 13/3*f1*x + 3/2*f1*x^2 - 1/6*f1*x^3 ```

integrate the polynomial over the full domain of the rule

```defint(P_of_x,'x',[0 5]) ```
```ans = 55/24*f4 + 5/24*f3 + 5/24*f2 + 55/24*f1 ```