.. _symbolic: ******************************** Using the Symbolic Math Toolbox ******************************** .. include:: replace.txt .. index:: symbolic math, math - analytic, math - symbolic The majority of our consideration relates to numerical algorithms. But here, we focus on analytic solutions to mathematical problems. We seek analytical solutions in mathematics courses, and we typically solve them with pencil and paper. Computers are, of course, more suited to numerical algorithms. But |M| does have the ability to find analytic solutions when they are available. When a computer solves a mathematics problem analytically, we say that the computer is operating on symbols in addition to numbers. With symbolic operations, the computer would simplify expressions such as :math:`6/3` to :math:`2`, but would leave numbers like :math:`\pi`, :math:`e`, and :math:`17/7` as irreducible numbers. The Symbolic MATH Toolbox (SMT) is required to access the features of |M| discussed here. The SMT is part of the student |M| bundle, but is an additional purchase with other |M| licenses. .. _ballup: Throwing a Ball Up =================== .. index:: subs, vpa As a quick tour of some of the features of the SMT, consider throwing a ball straight upward. The scene is depicted in :numref:`fig_ball`. The algebraic equations describing the force of gravity on a moving ball are readily available on the web or in physics and dynamics textbooks. But let us use the SMT to derive the equations from differential equations and apply them. In this example, we will use functions from the SMT with minimal comments about them. Then we will give more details about the SMT functions beginning in :ref:`symAlg`. .. _fig_ball: .. figure:: throwBall.png :align: center :width: 40% A ball is thrown upward from height :math:`y_0` with an initial velocity of :math:`V_0`. When the ball reaches its maximum height, :math:`y_{max}` at time :math:`t_{max}`, its velocity is zero, :math:`v(t_{max}) = 0`. The ball hits the ground at time :math:`t_{final}`, :math:`y(t_{final}) = 0`. The vertical position of the ball with respect to time is :math:`y(t)`. The velocity of the ball is the derivative of the position, .. math:: :label: eqn:velocity v(t) = \frac{dy(t)}{dt}. The derivative of the velocity is the acceleration, which in this case is a constant from gravity in the negative direction. The acceleration is the second derivative of the position. .. math:: :label: eqn:accel a(t) = \frac{dv(t)}{dt} = \frac{dy^2(t)}{dt^2} = g .. math:: g = -9.8 \; m/sec^2 \mbox{ or } g = -32 \; feet/sec^2 Our starting point, equation :eq:`eqn:accel`, is a second order differential equation. In addition to describing the differential equation, we need to describe the initial conditions of the problem. The initial height of the ball is given by :math:`y(0) = y_0`; and the initial velocity from the throw is :math:`v(0) = V_0`. Symbolic variables and functions are first created in |M| with the ``syms`` statement. Another way to create a symbolic variable is to pass a value or expression to the ``sym`` function, such as ``z = sym(5)``. The ``sym`` function returns a symbolic variable. The differential equations in equation :eq:`eqn:velocity` and equation :eq:`eqn:accel` are solved in the following |M| code. We will declare ``y(t)`` (position) as a symbolic function, and ``y0`` (initial position), ``V0`` (initial velocity), and ``g`` (gravity) as symbolic variables. Time, ``t``, is implicitly declared from ``y(t)``. |M| solves symbolic differential equations with the ``dsolve`` function. Note in the code below that the equal sign (``=``) is, as always, used for assignment statements, but the double equal sign (``==``) is used to declare equality relationships in equations. :: >> syms y(t) y0 V0 g >> diffeq = diff(y(t),t,2) == g diffeq = diff(y(t), t, 2) == g >> v(t) = diff(y(t),t) v(t) = diff(y(t), t) >> cond = [y(0) == y0, v(0) == V0] cond = [ y(0) == y0, subs(diff(y(t), t), t, 0) == V0] >> y(t) = dsolve(diffeq, cond) y(t) = (g*t^2)/2 + V0*t + y0 The ``diff`` function takes the derivative of the position to find the velocity. :: >> v(t) = diff(y(t),t) v(t) = V0 + g*t The maximum height of the ball occurs when the velocity is zero. The ``solve`` function finds solutions to algebraic equations. :: >> t_max = solve(v(t) == 0) t_max = -V0/g >> y_max = y(t_max) y_max = -V0^2/(2*g) + y0 The ball lands on the ground when ``y(t)`` is zero. Two answers are returned because ``y(t)`` is a quadratic equation. :: >> t_final = solve(y(t) == 0, t) t_final = -(V0 + (V0^2 - 2*g*y0)^(1/2))/g -(V0 - (V0^2 - 2*g*y0)^(1/2))/g We can give the variables numerical values and substitute (``subs``) the values into the equations, and see a decimal version of the number with ``vpa``. :: % From 1 meter off groud, throw up at 10 m/s >> y0 = 1; V0 = 10; g = -9.8; >> subs(t_max) % time to max height ans = 50/49 >> vpa(subs(t_max)) ans = 1.0204081632653061224489795918367 % seconds >> subs(y_max) ans = 299/49 >> vpa(subs(y_max)) ans = 6.1020408163265306122448979591837 % meters >> subs(t_final(1)) % time until the ball hits the ground ans = (5^(1/2)*598^(1/2))/49 + 50/49 >> vpa(subs(t_final(1))) ans = 2.1363447440401890728050707611418 % seconds .. _symAlg: Symbolic Algebra ==================== Collect ------- .. index:: collect The ``collect`` function combines variables of the same power of the specified variable in an expression. :: >> syms x y >> f = (2*x^2 + 3*y) * (x^3 + x^2 + x*y^2 + 2*y); >> collect(f,x) ans = 2*x^5 + 2*x^4 + (2*y^2 + 3*y)*x^3 + 7*y*x^2 + 3*y^3*x + 6*y^2 >> collect(f,y) ans = 3*x*y^3 + (2*x^3 + 6)*y^2 + (3*x^3 + 7*x^2)*y + 2*x^2*(x^3 + x^2) Factor ------- .. index:: factor The ``factor`` function finds the set of irreducible polynomials whose product form the higher order polynomial argument to ``factor``. The factoring operation is useful for finding the roots of the polynomial. :: >> g = x^4 + 3*x^3 - 19*x^2 - 27*x + 90; >> factor(g) ans = [ x - 2, x - 3, x + 5, x + 3] Expand ------- .. index:: expand The ``expand`` function multiplies polynomial products to get a single polynomial. :: >> f = (3*x^2 + 2*x)*(2*x - 5); >> expand(f) ans = 6*x^3 - 11*x^2 - 10*x >> g g = (x - 2)*(x - 3)*(x + 3)*(x + 5) >> expand(g) ans = x^4 + 3*x^3 - 19*x^2 - 27*x + 90 Simplify --------- .. index:: simplify The ``simplify`` function generates a simpler form of an expression. :: >> g g = x^4 + 3*x^3 - 19*x^2 - 27*x + 90 >> h = expand((x-2)*(x+5)) h = x^2 + 3*x - 10 >> k = g/h k = (x^4 + 3*x^3 - 19*x^2 - 27*x + 90)/(x^2 + 3*x - 10) >> simplify(k) ans = x^2 - 9 >> f = cos(x)*cos(y) + sin(x)*sin(y); >> simplify(f) ans = cos(x - y) >> f = cos(x)^2 + sin(x)^2; >> simplify(f) ans = 1 .. _solve: Solve ------- .. index:: solve .. index:: roots of a function, function roots, roots .. index:: roots of a polynomial, polynomial roots Finding the roots of equations, i.e., :math:`x`, where :math:`f(x) = 0`, is addressed here and also in the context of two other topics. * To find the numerical roots of a polynomial using a linear algebra algorithm, see the ``roots`` function described in :ref:`eigroots`. * To use numerical methods to find the roots of any real equation, including non-polynomials, see the ``fzero`` function described in :ref:`fzero`. The SMT function finds the roots of any function. The ``solve(s,x)`` function returns the set of symbolic values of ``x`` that satisfy the expression ``s``. The default is to solve for ``s == 0``, but an equality statement may be used as part of the expression. Be certain to use an equality statement (``==``), not an assignment statement (``=``). :: >> g(x) = x^2 - 1; >> solve(g(x),x) ans = -1 1 >> solve(g(x) == 4,x) ans = 5^(1/2) -5^(1/2) >> h(x) = x + 2; >> solve(g(x) == h(x), x) ans = 1/2 - 13^(1/2)/2 13^(1/2)/2 + 1/2 >> syms t y(t) >> y(t) = cos(pi*t) % non-polynomial y(t) = cos(pi*t) >> solve(y(t), t) ans = 1/2 >> y(t) = 1 + t^2; >> solve(y(t), t) ans = % even complex roots -1i 1i Subs ------- .. index:: subs The ``subs(s)`` function returns a copy of ``s`` after replacing the symbolic variables in ``s`` with their values from the |M| workspace. The ``subs`` function was demonstrated in :ref:`ballup`. :: >> syms g(x) a b x >> g(x) = a*x^2 + b*x; >> a = 2; b = 6; >> subs(g) ans(x) = 2*x^2 + 6*x >> subs(g(3)) ans = 36 >> syms c >> h = g/c h(x) = (a*x^2 + b*x)/c >> c = 3; >> subs(h(3)) ans = 12 Vpa --- .. index:: vpa The Symbolic Math Toolbox tries to keep variables as exact symbolic values rather than as decimal approximations. The ``vpa`` function converts symbolic expressions to decimal values. The ``vpa`` function gives decimal answers. The name ``vpa`` is an acronym for variable--precision floating-point arithmetic. As the name implies, the ``vpa(x, d)`` function evaluates ``x`` to at least ``d`` digits. The default number of digits is 32. The maximum number of digits is :math:`2^{29}`. :: >> g(x) = x^2 - 1; >> h = solve(g(x) == 4, x) h = 5^(1/2) -5^(1/2) >> vpa(h) ans = 2.2360679774997896964091736687313 -2.2360679774997896964091736687313 Symbolic Calculus ==================== .. _symDerivatve: Symbolic Derivatives -------------------- .. index:: derivatives - symbolic, diff The ``diff(s, var, n)`` function takes the symbolic derivative of an expression with respect to a specified variable. An optional third argument specifies the order of the derivative. As with many SMT functions, the default variable is ``x``. The variable is needed when the equation has more than one variable, but it never hurts to specify it. :: >> syms x >> diff(x^2 + 3*x + 2,x) ans = 2*x + 3 >> diff(x^2 + 3*x + 2,x,2) % 2nd derivative ans = 2 >> diff(x*sin(3*x)) ans = sin(3*x) + 3*x*cos(3*x) >> diff(exp(-x/2)) ans = -exp(-x/2)/2 Symbolic Integration --------------------- .. index:: integration - symbolic, int The ``int`` function performs integration of an expression. The variable specification is optional if the expression has only one symbolic variable. Thus valid forms for indefinite integrals are ``int(s)`` and ``int(s,var)``. .. math:: \int x^2 + 3\,x + 2 \, dx = \frac{1}{3} x^3 + \frac{3}{2} x^2 + 2\,x :: >> syms x t >> int(x^2 + 3*x + 2,x) ans = (x*(2*x^2 + 9*x + 12))/6 >> expand(int(x^2 + 3*x + 2,x)) ans = x^3/3 + (3*x^2)/2 + 2*x For definite integrals, we add the limits of integration to the arguments: ``int(s,var,a,b)``. .. math:: \int_0^5 x^2 + 3\,x + 2 \, dx = \frac{535}{6} :: >> int(x^2 + 3*x + 2,x,0,5) % definite integral ans = 535/6 .. math:: \int_0^t x^2 + 3\,x + 2 \, dx = \frac{1}{3} t^3 + \frac{3}{2} t^2 + 2\,t :: >> expand(int(x^2 + 3*x + 2,x,0,t)) ans = t^3/3 + (3*t^2)/2 + 2*t .. math:: \int_0^t e^{-x/2}\, dx = 2 - 2\,e^{-t/2} :: >> int(exp(-x/2),x,0,t) ans = 2 - 2*exp(-t/2) .. math:: \int_0^\infty e^{-x/2}\, dx = 2 :: >> int(exp(-x/2),x,0,Inf) % definite to infinity ans = 2 When the symbolic math toolbox is not able to find an analytic solution, it may return a result that makes use of numerical integration methods. We also encountered the ``erf`` and ``erfc`` functions in :ref:`continuous-distributions`. :: >> p = int(exp(-x^2), -Inf, 2) p = (pi^(1/2)*(erf(2) + 1))/2 Symbolic Limits --------------------- .. index:: limit - symbolic The ``limit(s, h, k)`` function returns the limit of an expression, ``s`` as the variable ``h`` approaches ``k``. .. math:: \lim_{x \rightarrow 0} \frac{\sin(x)}{x} = 1 .. math:: \lim_{x \rightarrow \infty} x\,e^{-x} = 0 :: >> limit(sin(x)/x,x,0) ans = 1 >> limit(x*exp(-x),x,Inf) ans = 0 Another optional argument of ``'left'`` or ``'right'`` is useful when the expression has a discontinuity at :math:`h = k`. The expression :math:`\frac{\abs{x - 3}}{x -3}` has a discontinuity at :math:`x = 3`. :: >> limit(abs(x-3)/(x-3),x,3,'right') ans = 1 >> limit(abs(x-3)/(x-3),x,3,'left') ans = -1 The following computes :math:`\frac{dy}{dx}` for :math:`y(x) = x^2` using a limit. Recall the definition of a derivative: .. math:: \frac{dy}{dx} = \lim_{h \rightarrow 0} \frac{y(x + h) - y(x)}{h}. :: >> syms h s x >> s = ((x + h)^2 - x^2)/h; >> limit(s,h,0) ans = 2*x .. _symbODE: Symbolic Differential Equations ================================ .. index:: differential equations - symbolic, dsolve, fplot Although most students in this course have not yet studied the analytic math techniques for solving differential equations, we use the solutions to differential equations in all areas of physics and engineering. Differential equations are just expressions where one or more terms of the expression is the derivative of another variable. Any system that moves or changes in either time or space is defined by a differential equation. The ``dsolve`` function returns the solution to differential equations. Note that previous versions of |M| used the special variables ``Dy`` and ``D2y`` for the first and second derivatives of :math:`y`. That notation is now deprecated. Instead use the ``diff`` function to define derivatives in the expression. As an example, consider the simple differential equation in :eq:`eqn-diff-exp`, which is solved by the following |M| code. .. math:: :label: eqn-diff-exp \begin{array}{l} \frac{dy(t)}{dt} = -a\,y(t) \\ \\ y(t) = C\,e^{-a\,t} \end{array} :: >> syms t y(t) a >> y(t) = diff(y(t)) == -a*y(t); >> Y(t) = dsolve(y(t)) Y(t) = C4*exp(-a*t) The example above shows two important points about differential equations. First, the solution contains an exponential function of the irrational number :math:`e`. This is because :math:`e^{a\,t}` is the only function for which the derivative is a multiple of the original function. Specifically, .. math:: :label: eqn-exp \frac{d\,e^{a\,t}}{dt} = a\,e^{a\,t}. Thus, exponential functions of :math:`e` show up as the solution to any differential equation of growth or decay where there is a linear relationship between the rate of change and the function's value. Secondly, the solution to equation :eq:`eqn-exp` contains a constant, ``C4``. This is normal when an initial condition or boundary condition was not specified. With an initial condition, we can either solve for ``C4`` algebraically, or use the condition as an input to ``dsolve``. For example, let's say that :math:`Y(0) = 5`. From our solution, :math:`Y(0) = C4`, thus :math:`C4 = 5`. The following code uses the SMT to solve a first order differential equation like :eq:`eqn-exp` with exponential decay. The solution is plotted in :numref:`fig-decay-eqn`. Initial or boundary conditions may also be passed to ``dsolve`` as an extra argument. When there is more than one condition, as is common for second or higher order differential equations, the conditions are passed as an array. See the :ref:`ballup` example for a second order differential equation example with two initial conditions. :: >> syms y(t) a >> eqn = diff(y(t)) == -a*y(t); >> cond = y(0) == 5; >> Y(t) = dsolve(eqn,cond) Y(t) = 5*exp(-a*t) .. note:: As with function handles, the results computed by the SMT that are functions may be plotted with the ``fplot`` function. .. index:: fplot :: >> a = 0.005; >> y(t) = subs(Y(t)) y(t) = 5*exp(-t/200) >> fplot(y(t), [0 1000]) .. _fig-decay-eqn: .. figure:: decayEqn.png :align: center :width: 60% The exponential decay solution to a differential equation of the same form as equation :eq:`eqn-exp`.