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
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)
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
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
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
% 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
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)
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)
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