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