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 (needed for Octave) pkg load symbolic syms initpython sympref display unicode
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.
disp(' ') % 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 disp('Integrate from 0 to 2 using "quad":') qNumInt02=quad(qNum,0,2) disp('We got a floating point number.') % find a root disp('Find one of the two roots using "fzero":') qNumRoot=fzero(qNum,1) disp('We got a floating point number.') % tell Matlab that, from now on, x is a symbolic variable disp(' ') disp('Let''s do some *symbolic math* now:') syms x % the example quadratic as a symbolic function disp('Create a symbolic quadratic:') qSym=-4*x^2+3*x+12 % we can integrate it between, say, 0 and 2 disp('Integrate symbolically from 0 to 2 using "int":') qSymInt02=int(qSym,x,[0 2]) disp('We got an *exact* number.') % but we can also find the anti-derivative disp('Find the symbolical anti-derivative using "int":') qSymInt=int(qSym,x) disp('We got a symbolic function, the anti-derivative.') disp('This is already in expanded form in Octave:') expand(qSymInt) % and we can find the derivative disp('Differentiate symbolically using "diff":') qSymDiff=diff(qSym,x) disp('We got a symbolic function, the derivative.') % we can find both roots exactly (note ==, not =) disp('Find both roots exactly using "solve" (note ==):') qSymRoots=solve(qSym == 0,x) disp('Display flat:') 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 disp('Try to factor the quadratic:') qSymFactors=factor(qSym) % Octave help on this is under @sym/ not sym/ like Matlab disp('Oops! Try "help @sym/factor", not "help factor".') disp('Note: in Octave "@sym/...", in Matlab "sym/..."') %help @sym/factor disp('Oh well. Octave is free software.') % we can easily factor a quadratic with rational roots disp('Try to factor a quadratic with rational roots:') qSymRat=-4*x^2+3*x+27 qSymRatFactors=factor(qSymRat)
qNum =
@(x) -4 * x .^ 2 + 3 * x + 12
Integrate from 0 to 2 using "quad":
qNumInt02 = 19.333
We got a floating point number.
Find one of the two roots using "fzero":
qNumRoot = -1.3972
We got a floating point number.
Let's do some *symbolic math* now:
Create a symbolic quadratic:
qSym = (sym)
2
- 4⋅x + 3⋅x + 12
Integrate symbolically from 0 to 2 using "int":
qSymInt02 = (sym) 58/3
We got an *exact* number.
Find the symbolical anti-derivative using "int":
qSymInt = (sym)
3 2
4⋅x 3⋅x
- ──── + ──── + 12⋅x
3 2
We got a symbolic function, the anti-derivative.
This is already in expanded form in Octave:
ans = (sym)
3 2
4⋅x 3⋅x
- ──── + ──── + 12⋅x
3 2
Differentiate symbolically using "diff":
qSymDiff = (sym) -8⋅x + 3
We got a symbolic function, the derivative.
Find both roots exactly using "solve" (note ==):
qSymRoots = (sym 2×1 matrix)
⎡ _____ ⎤
⎢ 3 ╲╱ 201 ⎥
⎢ ─ + ─────── ⎥
⎢ 8 8 ⎥
⎢ ⎥
⎢ _____ ⎥
⎢ ╲╱ 201 3⎥
⎢- ─────── + ─⎥
⎣ 8 8⎦
Display flat:
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
Try to factor the quadratic:
qSymFactors = (sym)
2
- 4⋅x + 3⋅x + 12
Oops! Try "help @sym/factor", not "help factor".
Note: in Octave "@sym/...", in Matlab "sym/..."
Oh well. Octave is free software.
Try to factor a quadratic with rational roots:
qSymRat = (sym)
2
- 4⋅x + 3⋅x + 27
qSymRatFactors = (sym) -(x - 3)⋅(4⋅x + 9)
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.6491... and gives that to the symbolic simplify
function. Of course simplify cannot make any sense
out of 0.6491.... 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.
disp(' ') % not so easy to see that 17 is a common factor disp('simplify(629/969) is *wrong*!') %ratioSimplified=simplify(629/969) disp('simplify(sym(''629/969'') is right:') ratioSimplified=simplify(sym('629/969'))
simplify(629/969) is *wrong*!
simplify(sym('629/969') is right:
ratioSimplified = (sym)
37
──
57
Function vpa, ("variable precision arithmetic"), will
give you numbers to arbitrarily high accuracy.
But watch it again. Let's try the result pi^2/6 of the infinite sum we talked about in the previous lesson.
disp(' ') % Octave 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
The Symbolic Toolbox can solve quite a lot of equations exactly.
If it cannot, it will drop back to a numerical solution.
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
disp(' ') % 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 disp('Find the roots of this quadratic using "solve"') 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 disp('Test it for -4x^2 + 3x +12 = 0:') qGenRootsTest=subs(qGenRoots,{a b c},{-4 3 12}) % we can create a normal Matlab function for the roots disp('Create a *numerical* function for the roots:') qGenRootsFun=matlabFunction(qGenRoots) % test it disp('Test that for -4x^2 + 3x +12 = 0:') qGenRootsFunTest=qGenRootsFun(-4,3,12)
qGen = (sym)
2
a⋅x + b⋅x + c
Find the roots of this quadratic using "solve"
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)
Test it for -4x^2 + 3x +12 = 0:
qGenRootsTest = (sym 2×1 matrix)
⎡ _____ ⎤
⎢ ╲╱ 201 3⎥
⎢- ─────── + ─⎥
⎢ 8 8⎥
⎢ ⎥
⎢ _____ ⎥
⎢ 3 ╲╱ 201 ⎥
⎢ ─ + ─────── ⎥
⎣ 8 8 ⎦
Create a *numerical* function for the roots:
qGenRootsFun =
@(a, b, c) [(-b + sqrt(-4 .* a .* c + b .^ 2)) ./ (2 .* a); -(b + sqrt(-4 .* a .* c + b .^ 2)) ./ (2 .* a)]
Test that for -4x^2 + 3x +12 = 0:
qGenRootsFunTest =
-1.3972
2.1472
Some other equations package Symbolic manages to solve.
disp(' ') % make sure x and y are symbols syms x y % exponentials can be negative for complex arguments disp('Solve e^x = -1:') expSol=solve(exp(x) == -1,x) % the next expression to solve, but we will write it out disp(' ') disp('Let''s try to solve a x y - b x - 1 = 0 now:') qxySym=a*x*y - b*x - 1 % a x y - b x - 1 = 0 can be solved for either x or y disp('Solve for x:') xSol=solve(a*x*y - b*x - 1 == 0, x) disp('Solve for y instead:') ySol=solve(a*x*y - b*x - 1 == 0, y) % sometimes rewriting the equation can help disp('Octave does not have the "collect" function.') %xCoefs=collect(a*x*y-b*x-1,x) % solve can solve systems of equations disp(' ') disp('Solve two linear equations in two unknowns:') [xSol,ySol] = solve(x+y == 7, x-y == 1,x,y) % if solve fails to do so, Octave fails disp(' ') disp('Try to solve two non-linear linear equations:') solve(x^2+cos(y) == 7, cosh(x)-y == 1,x,y) disp('Failed.')
Solve e^x = -1:
expSol = (sym) ⅈ⋅π
Let's try to solve a x y - b x - 1 = 0 now:
qxySym = (sym) a⋅x⋅y - b⋅x - 1
Solve for x:
xSol = (sym)
1
───────
a⋅y - b
Solve for y instead:
ySol = (sym)
b⋅x + 1
───────
a⋅x
Octave does not have the "collect" function.
Solve two linear equations in two unknowns:
xSol = (sym) 4
ySol = (sym) 3
Try to solve two non-linear linear equations:
ans = {}(0x0)
Failed.
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.
disp(' ') % make sure x is still symbolic syms x % example symbolic ratio disp('Create an 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 disp('Let''s look at the factors first:') ratSymFactors=factor(ratSym) disp('Already multiplied out.') %prod(ratSymFactors) % get the partial fractions disp('Now get the partial fractions') ratPartFrac=partfrac(ratSym)
Create an example symbolic ratio:
ratSym = (sym)
2
2⋅x - 3⋅x + 1
────────────────────
3 2
x + 2⋅x - 9⋅x - 18
Let's look at the factors first:
ratSymFactors = (sym)
(x - 1)⋅(2⋅x - 1)
───────────────────────
(x - 3)⋅(x + 2)⋅(x + 3)
Already multiplied out.
Now get the partial fractions
ratPartFrac = (sym)
14 3 1
───────── - ───── + ─────────
3⋅(x + 3) x + 2 3⋅(x - 3)
Below are some more example manipulations
disp(' ') % make sure x, etc, are still symbolic syms x a b c % dealing with integrals you cannot do disp(' ') disp('Try two hard integrals from a basic table book:') disp('Try 1/(x*(a*x^2+b*x+c)^(1/2)):') intHard1=int(sym('1/(x*(a*x^2+b*x+c)^(1/2))')) disp('Try 1/(x*(a*x^2+b*x+c)^(3/2)):') intHard2=int(sym('1/(x*(a*x^2+b*x+c)^(3/2))')) disp('But the solution is in a basic table book!') disp('Well, Octave is free (and so is Sympy).') % sometimes int is useful to identify a function disp(' ') disp('The integral of exp(-x^2):') intExpMinusxSqr=int(exp(-x^2),x) disp('The integral of sin(x)/x:') intSinxOverx=int(sin(x)/x,x) % writing the Taylor series of, say, sinint disp(' ') disp('Try to learn about the Taylor series of Si:') sinxOverxTaylor=taylor(sin(x)/x) sinxOverxTaylor=taylor(sin(x)/x,'Order',10) sinintTaylor=int(sinxOverxTaylor) % examine the logic of the terms disp('Examine t(3)/t(1):') t3ot1=simplify(-(x^3/18)/x) factors18=factor(18) disp('Examine t(5)/t(3):') t5ot3=simplify(-(x^5/600)/(x^3/18)) factors100=factor(100) disp('Examine t(7)/t(5):') t7ot5=simplify(-(x^7/35280)/(x^5/600)) factors294=factor(294) disp('Examine t(9)/t(7):') t9ot7=simplify(-(x^9/3265920)/(x^7/35280)) factors648=factor(648) % funtool has not yet been implemented in Octave
Try two hard integrals from a basic table book:
Try 1/(x*(a*x^2+b*x+c)^(1/2)):
intHard1 = (sym)
⌠
⎮ 1
⎮ ───────────────────── dx
⎮ ________________
⎮ ╱ 2
⎮ x⋅╲╱ a⋅x + b⋅x + c
⌡
Try 1/(x*(a*x^2+b*x+c)^(3/2)):
intHard2 = (sym)
⌠
⎮ 1
⎮ ───────────────────── dx
⎮ 3/2
⎮ ⎛ 2 ⎞
⎮ x⋅⎝a⋅x + b⋅x + c⎠
⌡
But the solution is in a basic table book!
Well, Octave is free (and so is Sympy).
The integral of exp(-x^2):
intExpMinusxSqr = (sym)
___
╲╱ π ⋅erf(x)
────────────
2
The integral of sin(x)/x:
intSinxOverx = (sym) Si(x)
Try to learn about the Taylor series of Si:
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):
t3ot1 = (sym)
2
-x
────
18
factors18 =
2 3 3
Examine t(5)/t(3):
t5ot3 = (sym)
2
-3⋅x
──────
100
factors100 =
2 2 5 5
Examine t(7)/t(5):
t7ot5 = (sym)
2
-5⋅x
──────
294
factors294 =
2 3 7 7
Examine t(9)/t(7):
t9ot7 = (sym)
2
-7⋅x
──────
648
factors648 =
2 2 2 3 3 3 3