7 SYMBOLIC COMPUTATIONS

Contents

Initialization

Note: Either publish to html or use

sympref display ascii

not unicode. Latex does not have various Unicode symbols implemented.

% reduce needless whitespace
format compact
% reduce irritations
more off
% start a diary
%diary lectureN.txt

% load the symbolic package
pkg load symbolic
syms initpython
sympref display unicode

INTRODUCTION

Normally Matlab generates numbers, not formulae. If you ask Matlab to find the roots where a function is zero, maybe using fzero, it gives you numbers. If you ask Matlab to integrate a function from a lower limit to an upper limit using integral or quad, it gives you a number. Etcetera,

But sometimes you want a formula instead of a number. For example, you might want the derivative or antiderivative of a function. Either one is a formula, not a number. Also, sometimes you want to see an exact number, and non-integer numbers in Matlab have round-off errors.

What you need for such purposes is a symbolic computations program. Such a program is available inside Matlab as the "Symbolic Math Toolbox". Normally MathWorks charges separately for this package. However, it is included in the "Student Edition", and within Matlab on the COE computers.

In this section we will illustrate how normal Matlab computations and symbolic ones differ. We will look at a simple function, a quadratic in fact.

One thing to remember: Be sure to inform Matlab with the syms or sym command when variables and/or numbers are intended to be symbolic. Normal variables are names of storage locations with a number in it. But a symbolic variable does not store a number; at all times it can stand for any number. So a normal variable x is very different from a symbolic variable x, and Matlab must know which of the two x is.

% the example quadratic as a normal Matlab function
qNum = @(x) -4*x.^2+3*x+12

% we can integrate it between, say, 0 and 2
qNumInt02=quad(qNum,0,2)
disp('We got a floating point number.')

% find a root
qNumRoot=fzero(qNum,1)
disp('We got a floating point number.')

% tell Matlab that, from now on, x is a symbolic variable
syms x

% the example quadratic now as a symbolic function
qSym=-4*x^2+3*x+12

% we can integrate it between, say, 0 and 2
qSymInt02=int(qSym,x,[0 2])
disp('We got an exact number.')

% but we can also find the anti-derivative
qSymInt=int(qSym,x)
disp('We got a function, the anti-derivative.')
disp('This is already in expanded form in Octave:')
expand(qSymInt)

% and we can find the derivative
qSymDiff=diff(qSym,x)
disp('We got a function, the derivative.')

% we can find both roots (note ==, not =)
qSymRoots=solve(qSym == 0,x)
sympref display flat
qSymRoots
sympref display unicode
disp('We got both roots as symbolic expressions.')
disp('To see the numbers, use "double":')
double(qSymRoots)

% we cannot factor this quadratic on octave
qSymFactors=factor(qSym)
% note that Octave help is under @sym/ not sym/
disp('Oops! Try "help @sym/factor", not "help factor".')
disp('Note: in Octave "@sym/...", in Matlab "sym/..."')
%help @sym/factor

% we can easily factor a quadratic with rational roots
qSymRat=-4*x^2+3*x+27
qSymRatFactors=factor(qSymRat)
qNum =
@(x) -4 * x .^ 2 + 3 * x + 12
qNumInt02 =  19.333
We got a floating point number.
qNumRoot = -1.3972
We got a floating point number.
qSym = (sym)
       2           
  - 4⋅x  + 3⋅x + 12
qSymInt02 = (sym) 58/3
We got an exact number.
qSymInt = (sym)
       3      2       
    4⋅x    3⋅x        
  - ──── + ──── + 12⋅x
     3      2         
We got a function, the anti-derivative.
This is already in expanded form in Octave:
ans = (sym)
       3      2       
    4⋅x    3⋅x        
  - ──── + ──── + 12⋅x
     3      2         
qSymDiff = (sym) -8⋅x + 3
We got a function, the derivative.
qSymRoots = (sym 2×1 matrix)
  ⎡       _____ ⎤
  ⎢ 3   ╲╱ 201  ⎥
  ⎢ ─ + ─────── ⎥
  ⎢ 8      8    ⎥
  ⎢             ⎥
  ⎢    _____    ⎥
  ⎢  ╲╱ 201    3⎥
  ⎢- ─────── + ─⎥
  ⎣     8      8⎦
qSymRoots = (sym) Matrix([[3/8 + sqrt(201)/8], [-sqrt(201)/8 + 3/8]])  (2x1 matrix)
We got both roots as symbolic expressions.
To see the numbers, use "double":
ans =
   2.1472
  -1.3972
qSymFactors = (sym)
       2           
  - 4⋅x  + 3⋅x + 12
Oops! Try "help @sym/factor", not "help factor".
Note: in Octave "@sym/...", in Matlab "sym/..."
qSymRat = (sym)
       2           
  - 4⋅x  + 3⋅x + 27
qSymRatFactors = (sym) -(x - 3)⋅(4⋅x + 9)

SIMPLIFYING ANSWERS

In classes like Analysis in Mechanical Engineering, you are required to simplify your answers. Symbolic math to the rescue!

Watch it, however. If you try to simplify, say, a numeric ratio like 629/969 as

 simplify(629/969)

then Matlab sees "629/969", evaluates that as 0.9491... and gives that to the symbolic simplify function. Of course simplify cannot make any sense out of 0.9491.... In fact, it will refuse to cooperate. What you need to do is tell Matlab that the entire "629/969" is to be treated as a symbolic expression, to be given to simplify "as is". You can do that with the sym function.

% not so easy to see that 17 is a common factor
%ratioSimplified=simplify(629/969)
ratioSimplified=simplify(sym('629/969'))
ratioSimplified = (sym)
  37
  ──
  57

VERY HIGH ACCURACY

Function vpa, ("variable precision arithmetic"), will give you numbers to arbitrarily high accuracy.

But watch it again.

% Matlab gives vpa a number equal to pi^2/6 to 16 digits:
disp('vpa(pi^2/6) is wrong:')
piSqOver6=vpa(pi^2/6,50)

% Matlab gives vpa the symbolic string pi^2/6:
disp('vpa(sym(''pi^2/6'')) is correct:')
piSqOver6=vpa(sym('pi^2/6'),50)
vpa(pi^2/6) is wrong:
piSqOver6 = (sym) 1.6449340668482264060656916626612655818462371826172
vpa(sym('pi^2/6')) is correct:
piSqOver6 = (sym) 1.6449340668482264364724151666460251892189499012068

SOLVING EQUATIONS

The Symbolic Toolbox can solve quite a lot of equations exactly.

If it cannot, it will drop back to a numerical solution.

Solving quadratic equations

Suppose you no longer remembered the solution to the quadratic equation

$$a x^2 + b x + c = 0$$

The symbolic toolbox can give it to you

% tell Matlab that the variables are symbolic
syms a b c x qGen

% define the generic quadratic polynomial
qGen=a*x^2+b*x+c

% find the symbolic solution
qGenRoots=solve(qGen == 0,x)
disp('Octave already used "pretty" on the roots:')
pretty(qGenRoots)
disp('but we can display them "flat":')
sympref display flat
qGenRoots
sympref display unicode

% see whether we get back our previous solution
qGenRootsTest=subs(qGenRoots,{a b c},{-4 3 12})

% we can create a normal Matlab function for the roots
qGenRootsFun=matlabFunction(qGenRoots)

% test it
qGenRootsFunTest=qGenRootsFun(-4,3,12)
qGen = (sym)
     2          
  a⋅x  + b⋅x + c
qGenRoots = (sym 2×1 matrix)
  ⎡         _____________  ⎤
  ⎢        ╱           2   ⎥
  ⎢ -b + ╲╱  -4⋅a⋅c + b    ⎥
  ⎢ ─────────────────────  ⎥
  ⎢          2⋅a           ⎥
  ⎢                        ⎥
  ⎢ ⎛       _____________⎞ ⎥
  ⎢ ⎜      ╱           2 ⎟ ⎥
  ⎢-⎝b + ╲╱  -4⋅a⋅c + b  ⎠ ⎥
  ⎢────────────────────────⎥
  ⎣          2⋅a           ⎦
Octave already used "pretty" on the roots:
  ⎡         _____________  ⎤
  ⎢        ╱           2   ⎥
  ⎢ -b + ╲╱  -4⋅a⋅c + b    ⎥
  ⎢ ─────────────────────  ⎥
  ⎢          2⋅a           ⎥
  ⎢                        ⎥
  ⎢ ⎛       _____________⎞ ⎥
  ⎢ ⎜      ╱           2 ⎟ ⎥
  ⎢-⎝b + ╲╱  -4⋅a⋅c + b  ⎠ ⎥
  ⎢────────────────────────⎥
  ⎣          2⋅a           ⎦
but we can display them "flat":
qGenRoots = (sym) Matrix([[(-b + sqrt(-4*a*c + b**2))/(2*a)], [-(b + sqrt(-4*a*c + b**2))/(2*a)]])  (2x1 matrix)
qGenRootsTest = (sym 2×1 matrix)
  ⎡    _____    ⎤
  ⎢  ╲╱ 201    3⎥
  ⎢- ─────── + ─⎥
  ⎢     8      8⎥
  ⎢             ⎥
  ⎢       _____ ⎥
  ⎢ 3   ╲╱ 201  ⎥
  ⎢ ─ + ─────── ⎥
  ⎣ 8      8    ⎦
qGenRootsFun =
@(a, b, c) [(-b + sqrt(-4 .* a .* c + b .^ 2)) ./ (2 .* a); -(b + sqrt(-4 .* a .* c + b .^ 2)) ./ (2 .* a)]
qGenRootsFunTest =
  -1.3972
   2.1472

Solving other equations

Some other equations Matlab manages to solve

% make sure x and y are symbols
syms x y

% exponentials can be negative for complex arguments
expSol=solve(exp(x) == -1,x)

% the next expression to solve, but we will write it out
qxySym=a*x*y - b*x - 1

% a x y - b x - 1 can be solved for either x or y
xSol=solve(a*x*y - b*x - 1 == 0, x)
ySol=solve(a*x*y - b*x - 1 == 0, y)

% sometimes rewriting the equation can help
%xCoefs=collect(5*x*y-9*x-1,x)
disp('Octave does not have the "collect" function.')

% solve can solve systems of equations
[xSol,ySol] = solve(x+y == 7, x-y == 1',x,y)

% if it cannot find an exact solution, Octave fails
solve(x^2+cos(y) == 7, cosh(x)-y == 1,x,y)
expSol = (sym) ⅈ⋅π
qxySym = (sym) a⋅x⋅y - b⋅x - 1
xSol = (sym)
     1   
  ───────
  a⋅y - b
ySol = (sym)
  b⋅x + 1
  ───────
    a⋅x  
Octave does not have the "collect" function.
xSol = (sym) 4
ySol = (sym) 3
ans = {}(0x0)

PARTIAL FRACTIONS

In analyzing the dynamics of controlled systems, you often encounter ratios of big polynomials. Then you want to take these ratios apart in simpler fractions. The reason is that these simpler fractions tell you many important things. For example, they tell you whether, if the system is disturbed, it will return to its normal position, versus, say, crash. And, if so, they will also tell you how fast the system will return to normal.

Symbolic function partfrac can do it.

% make sure x is still symbolic
syms x

% example symbolic ratio
ratSym=(2*x^2-3*x+1)/(x^3+2*x^2-9*x-18)

% you could first look at the factors
ratSymFactors=factor(ratSym)

% get the partial fractions
ratPartFrac=partfrac(ratSym)
ratSym = (sym)
        2             
     2⋅x  - 3⋅x + 1   
  ────────────────────
   3      2           
  x  + 2⋅x  - 9⋅x - 18
ratSymFactors = (sym)
     (x - 1)⋅(2⋅x - 1)   
  ───────────────────────
  (x - 3)⋅(x + 2)⋅(x + 3)
ratPartFrac = (sym)
      14        3         1    
  ───────── - ───── + ─────────
  3⋅(x + 3)   x + 2   3⋅(x - 3)

FUNCTION MANIPULATIONS

Below are some more example manipulations

% make sure x, etc, are still symbolic
syms x a b c

% dealing with integrals you cannot do
intHard1=int(sym('1/(x*(a*x^2+b*x+c)^(1/2))'))
%intHard2=int(sym('1/(x*(a*x^2+b*x+c)^(3/2))'))
disp('Well, Octave is free (and so is Sympy).')

% sometimes int is useful to identify a function
intUnk=int(sin(x)/x,x)

% writing the Taylor series of, say, sinint
sinxOverxTaylor=taylor(sin(x)/x)
sinxOverxTaylor=taylor(sin(x)/x,'Order',10)
sinintTaylor=int(sinxOverxTaylor)

% see the logic
disp('Examine t(3)/t(1):')
simplify(-(x^3/18)/x)
factor(18)
disp('Examine t(5)/t(3):')
simplify(-(x^5/600)/(x^3/18))
factor(100)
disp('Examine t(7)/t(5):')
simplify(-(x^7/35280)/(x^5/600))
factor(294)
disp('Examine t(9)/t(7):')
simplify(-(x^9/3265920)/(x^7/35280))
factor(648)

% funtool has not yet been implemented in Octave
intHard1 = (sym)
  ⌠                         
  ⎮           1             
  ⎮ ───────────────────── dx
  ⎮      ________________   
  ⎮     ╱    2              
  ⎮ x⋅╲╱  a⋅x  + b⋅x + c    
  ⌡                         
Well, Octave is free (and so is Sympy).
intUnk = (sym) Si(x)
sinxOverxTaylor = (sym)
    4    2    
   x    x     
  ─── - ── + 1
  120   6     
sinxOverxTaylor = (sym)
     8       6      4    2    
    x       x      x    x     
  ────── - ──── + ─── - ── + 1
  362880   5040   120   6     
sinintTaylor = (sym)
      9        7      5    3    
     x        x      x    x     
  ─────── - ───── + ─── - ── + x
  3265920   35280   600   18    
Examine t(3)/t(1):
ans = (sym)
    2 
  -x  
  ────
   18 
ans =
   2   3   3
Examine t(5)/t(3):
ans = (sym)
      2 
  -3⋅x  
  ──────
   100  
ans =
   2   2   5   5
Examine t(7)/t(5):
ans = (sym)
      2 
  -5⋅x  
  ──────
   294  
ans =
   2   3   7   7
Examine t(9)/t(7):
ans = (sym)
      2 
  -7⋅x  
  ──────
   648  
ans =
   2   2   2   3   3   3   3

End lesson 7