From 02869d18ff55bd35d4e389643dcdca2f06ecf39c Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Mon, 6 Feb 2023 15:46:24 +0100 Subject: [PATCH 1/8] Update SymPy chapter for 2023. - Improved wording in various places. - Added more hyperlinks to connect to mentioned topics. - Added an exercise to almost all sections. - Expanded some examples. --- sympy.rst | 342 ++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 268 insertions(+), 74 deletions(-) diff --git a/sympy.rst b/sympy.rst index 5095379d..3ebeee53 100644 --- a/sympy.rst +++ b/sympy.rst @@ -8,29 +8,43 @@ SymPy :jupyter-download-script:`sympy` or Jupyter notebook: :jupyter-download-notebook:`sympy`. -SymPy_ is an open source, collaboratively developed computer algebra system +Learning Objectives +=================== + +After completing this chapter readers will be able to: + +- Write mathematical expressions with symbols and functions using SymPy. +- Differentiate mathematical expressions using SymPy. +- Evaluate mathematical expressions using SymPy. +- Solve a linear system of equations with SymPy. + +Introduction +============ + +SymPy_ is an open source, collaboratively developed `computer algebra system`_ (CAS) written in Python. It will be used extensively for manipulating symbolic -expressions and equations. All of the mathematics needed to formulate the +expressions and equations. All of the mathematics needed to formulate the equations of motion of multibody systems can be done with pencil and paper, but -the bookkeeping becomes extremely tedious and error prone for systems with -even a small number of bodies. SymPy lets a computer handle the -tedious aspects, e.g. differentiation, solving linear systems of equations, and -reduces the errors one would encounter with pencil and paper. This chapter -introduces SymPy and the primary core SymPy features needed we will be using. +the bookkeeping becomes extremely tedious and error prone for systems with even +a small number of bodies. SymPy lets a computer handle the tedious aspects +(e.g. differentiation or solving linear systems of equations) and reduces the +errors one would encounter with pencil and paper. This chapter introduces SymPy +and the primary SymPy features needed we will be using. .. _SymPy: https://www.sympy.org +.. _computer algebra system: https://en.wikipedia.org/wiki/Computer_algebra_system Import and Setup ================ -I will consistently import SymPy as follows: +I will import SymPy as follows through this book: .. jupyter-execute:: import sympy as sm -Since SymPy works with mathematical symbols it's nice to view SymPy -objects in a format that is similar to the math in a textbook. Executing +Since SymPy works with mathematical symbols it's nice to view SymPy objects in +a format that is similar to the math in a textbook. Executing :external:py:func:`~sympy.interactive.printing.init_printing` at the beginning of your Jupyter Notebook will ensure that SymPy objects render as typeset mathematics. I use the ``use_latex='mathjax'`` argument here to disable math @@ -43,8 +57,9 @@ image generation. Symbols ======= -Symbols are created with the :external:py:func:`~sympy.core.symbol.symbols` -function. A symbol :math:`a` is created like so: +Mathematical symbols are created with the +:external:py:func:`~sympy.core.symbol.symbols` function. A symbol :math:`a` is +created like so: .. jupyter-execute:: @@ -62,8 +77,8 @@ recognizes common Greek symbols by their spelled-out name. .. jupyter-execute:: - b, t, omega = sm.symbols('b, t, omega') - b, t, omega + b, t, omega, Omega = sm.symbols('b, t, omega, Omega') + b, t, omega, Omega Note that the argument provided to ``symbols()`` does not need to match the Python variable name it is assigned to. Using more verbose Python variable @@ -76,6 +91,19 @@ are recognized too. pivot_angle, w2 = sm.symbols('alpha1, omega2') pivot_angle, w2 +.. admonition:: Exercise + + Review the SymPy documentation and create symbols :math:`q_1, q_2, \ldots, + q_{10}` with a very succint call to + :external:py:func:`~sympy.core.symbol.symbols`. + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + sm.symbols('q1:11') + Undefined Functions =================== @@ -102,14 +130,16 @@ Now you can create functions of one or more variables like so: f(t) -Due to SymPy's internal implementations, the type of this object is not defined -as expected: +.. warning:: -.. jupyter-execute:: + Due to SymPy's internal implementations, the type of a function with its + argument is not defined as expected: - type(f(t)) + .. jupyter-execute:: + + type(f(t)) -so be aware of that. + This can be confusing if you are checking types. The same ``UndefinedFunction`` can be used to create multivariate functions: @@ -117,6 +147,19 @@ The same ``UndefinedFunction`` can be used to create multivariate functions: f(a, b, omega, t) +.. admonition:: Exercise + + Create a function :math:`H(x, y, z)`. + + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + x, y, z = sm.symbols('x, y, z') + sm.Function('H')(x, y, z) + Symbolic Expressions ==================== @@ -136,8 +179,8 @@ An expression will have the type ``Add, Mul, or Pow``: type(expr1) -This is because SymPy stores expressions as a tree_. You can inspect this -internal representation by using the +This is because SymPy stores expressions behind the scenes as a tree_. You can +inspect this internal representation by using the :external:py:func:`~sympy.printing.repr.srepr` function: .. _tree: https://en.wikipedia.org/wiki/Tree_(graph_theory) @@ -193,7 +236,7 @@ the ``srepr()`` version can help clear up misunderstandings. See the .. _manipulation section: https://docs.sympy.org/latest/tutorial/manipulation.html -Undefined functions can also be used in expressions: +Undefined functions can also be used in expressions just like symbols: .. jupyter-execute:: @@ -239,6 +282,12 @@ expressions: sm.S(1)/2*a + Or you can ensure the symbol comes first in the division operation: + + .. jupyter-execute:: + + a/2 + Lastly, an expression of ``t``: .. jupyter-execute:: @@ -246,12 +295,31 @@ Lastly, an expression of ``t``: expr5 = t*sm.sin(omega*f(t)) + f(t)/sm.sqrt(t) expr5 +.. admonition:: Exercise + + Create an expression for the normal distribution function: + + .. math:: + + \frac{1}{\sqrt{2\pi\sigma}}e^{\frac{(x-\mu)^2}{2\sigma^2}} + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + x, s, m = sm.symbols('x, sigma, mu') + sm.exp((x-m)**2/2/s**2)/sm.sqrt(2*sm.pi*s) + + Notice that SymPy does some minor manipulation of the expression, but it is + equivalent to the form shown in the prompt. + Printing ======== -I introduced the ``srepr()`` form of SymPy expressions above and mentioned -that expressions can have different representations. For the following -``srepr()`` form: +I introduced the ``srepr()`` form of SymPy expressions above and mentioned that +expressions can have different representations. For the following ``srepr()`` +form: .. jupyter-execute:: @@ -279,8 +347,9 @@ to provide a form that more closely resembles typeset math: sm.pprint(expr3) Lastly, the following lines show how SymPy expressions can be represented as -LaTeX code. The double backslashes are present because double backslashes -represent the escape character in Python strings. +LaTeX code using :external:py:func:`sympy.printing.latex.latex`. The double +backslashes are present because double backslashes represent the escape +character in Python strings. .. jupyter-execute:: @@ -297,6 +366,25 @@ represent the escape character in Python strings. to the screen make take a long time and fill your entire notebook with an unreadable mess. +.. admonition:: Exercise + + Print the normal distribution expression + + .. math:: + + \frac{1}{\sqrt{2\pi\sigma}}e^{\frac{(x-\mu)^2}{2\sigma^2}} + + as a LaTeX string inside an equation environment. + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + x, s, m = sm.symbols('x, sigma, mu') + print(sm.latex(sm.exp((x-m)**2/2/s**2)/sm.sqrt(2*sm.pi*s), + mode='equation')) + Differentiating =============== @@ -330,12 +418,21 @@ It can be differentiated, for example, with respect to :math:`b`: expr3.diff(b) You can also calculate partial derivatives with respect to successive -variables, as in the following: +variables. If you want to first differentiate with respect to :math:`b` and +then with respect to :math:`t` as in the following operation: .. math:: \frac{\partial^2 h(a, \omega, t, b)}{\partial t \partial b} +where: + +.. math:: + + h(a, \omega, t, b) = \displaystyle a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}} + +then you can use successive arguments to ``.diff()``: + .. jupyter-execute:: expr3.diff(b, t) @@ -367,18 +464,37 @@ function`_. sm.Abs(h(t)).diff(t) Sometimes you may need to add assumptions to variables, but in general it - will not be necessary. + will not be necessary. Read more about `assumptions in SymPy's guide + `_. -Lastly, a typical type of derivative you may encounter: +.. admonition:: Exercise -.. jupyter-execute:: + Differentiate ``expr5`` above using this operator: - expr5 + .. math:: -.. jupyter-execute:: + \frac{\partial^2}{\partial \omega \partial t} + +.. admonition:: Solution + :class: dropdown + + First show ``expr5``: + + .. jupyter-execute:: + + expr5 - expr5.diff(t) + The twice partial derivative is: + + .. jupyter-execute:: + + expr5.diff(t, omega) + + or you can chain ``.diff()`` calls: + + .. jupyter-execute:: + expr5.diff(t).diff(omega) Evaluating Symbolic Expressions =============================== @@ -406,19 +522,19 @@ the number of decimal places and the substitution dictionary to evaluate: .. jupyter-execute:: - expr3.evalf(n=31, subs=repl) + expr3.evalf(n=10, subs=repl) .. jupyter-execute:: - type(expr3.evalf(n=31, subs=repl)) + type(expr3.evalf(n=10, subs=repl)) -Note that this is a SymPy :external:py:class:`~sympy.core.numbers.Float` object, which is a special object that can -have an arbitrary number of decimal places, for example here is the expression -evaluated to 300 decimal places: +Note that this is a SymPy :external:py:class:`~sympy.core.numbers.Float` +object, which is a special object that can have an arbitrary number of decimal +places, for example here is the expression evaluated to 80 decimal places: .. jupyter-execute:: - expr3.evalf(n=300, subs=repl) + expr3.evalf(n=80, subs=repl) To convert this to Python floating point number, use ``float()``: @@ -430,40 +546,95 @@ To convert this to Python floating point number, use ``float()``: type(float(expr3.evalf(n=300, subs=repl))) -This value is a machine precision floating point value and can be used with +This value is a `machine precision`_ floating point value and can be used with standard Python functions that operate on floating point numbers. -To obtain machine precision floating point numbers directly, it is better to -use the :external:py:func:`~sympy.utilities.lambdify.lambdify` function to convert the expression to a Python -function: +.. _machine precision: https://en.wikipedia.org/wiki/Machine_epsilon + +To obtain machine precision floating point numbers directly and with more +flexibility, it is better to use the +:external:py:func:`~sympy.utilities.lambdify.lambdify` function to convert the +expression to a Python function. When using ``lambdify()``, all symbols and +functions should be converted to numbers so first identify what symbols and +functions make up the expression. + +.. jupyter-execute:: + + expr3 + +:math:`\omega, a, f(t)`, and :math:`b` are all present in the expression. The +first argument of ``lambdify()`` should be a sequence of all these symbols and +functions and the second argument should be the expression. .. jupyter-execute:: eval_expr3 = sm.lambdify((omega, a, f(t), b), expr3) +``lambdify()`` generates a Python function and, in this case, we store that +function in the variable ``eval_expr3``. You can see what the inputs and +outputs of the function are with ``help()``: + .. jupyter-execute:: help(eval_expr3) -Now you have a function that operates on and returns floating point values: +This function operates on and returns floating point values, for example: .. jupyter-execute:: eval_expr3(3.14/4, 2, -12, 25) +The type of lambdify's return values will be NumPy_ floats. + .. jupyter-execute:: type(eval_expr3(3.14/4, 2, -12, 25)) -This distinction between SymPy ``Float`` objects and regular Python and NumPy -``float`` objects is important. In this case, the Python float and the NumPy -float are equivalent. The later will compute much faster because arbitrary -precision is not required. +.. _NumPy: https://www.numpy.org + +These floats are interoperable with Python floats for single values (unlike +SymPy Floats) but also support arrays of floats. For example: + +.. jupyter-execute:: + + eval_expr3(3.14/4, 2, -12, [25, 26, 27]) + +.. jupyter-execute:: + + type(eval_expr3(3.14/4, 2, -12, [25, 26, 27])) + +More on NumPy arrays of floats will be introduced in a later chapter. + +.. warning:: Python and NumPy floats can be mixed, but avoid mixing SymPy + Floats with either. .. note:: - In this course, you will almost always want to convert SymPy expressions - to machine precision floating point numbers, so use ``lambdify()``. + This distinction between SymPy ``Float`` objects and regular Python and + NumPy ``float`` objects is important. In this case, the Python float and the + NumPy float are equivalent. The later will compute much faster because + arbitrary precision is not required. In this course, you will almost always + want to convert SymPy expressions to machine precision floating point + numbers, so use ``lambdify()``. + +.. admonition:: Exercise + + Create a symbolic expression representing `Newton's Law of Universal + Gravitation + `_. Use + ``lambdify()`` to evaluate the expression for two mass of 5.972E24 kg and 80 + kg at a distance of 6371 km apart to find the gravitational force in + Newtons. + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + G, m1, m2, r = sm.symbols('G, m1, m2, r') + F = G*m1*m2/r**2 + eval_F = sm.lambdify((G, m1, m2, r), F) + eval_F(6.67430E-11, 5.972E24, 80, 6371E3) Matrices ======== @@ -471,7 +642,8 @@ Matrices SymPy supports matrices of expressions and linear algebra. Many of the operations needed in multibody dynamics are more succinctly formulated with matrices and linear algebra. Matrices can be created by passing nested lists to -the :external:py:class:`Matrix() ` object. For example: +the :external:py:class:`Matrix() ` +object. For example: .. jupyter-execute:: @@ -548,7 +720,11 @@ Differentiation operates on each element of the matrix: mat3.diff(t) -The Jacobian_ matrix of vector (column matrix) can be formed with the +If you have column vectors :math:`\bar{v}` and :math:`\bar{u}`, the +:math:`(i,j)` entries of the Jacobian of :math:`\bar{v}` with respect to the +entries in vector :math:`\bar{u}` are found with :math:`\mathbf{J}_{ij} = +\frac{\partial v_i}{\partial u_j}`. The Jacobian_ matrix of vector (column +matrix) can be formed with the :external:py:meth:`~sympy.matrices.matrices.MatrixCalculus.jacobian` method. This calculates the partial derivatives of each element in the vector with respect to a vector (or sequence) of variables. @@ -564,21 +740,25 @@ respect to a vector (or sequence) of variables. .. _Jacobian: https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant +.. todo:: Add an exercise here. Solving Linear Systems ====================== -You'll need to solve linear systems of equations often in this course. SymPy -offers a number of ways to do this, but the best way to do so if you know a set -of equations are linear in specific variables is the method described below. -First, you should know you have equations of this form: +You'll need to solve `linear systems of equations`_ often in this course. SymPy +offers `a number of ways to do this`_, but the best way to do so if you know a +set of equations are linear in specific variables is the method described +below. First, you should know you have equations of this form: .. math:: - a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n + b_1 = 0 \\ - a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n + b_2 = 0 \\ - \ldots \\ - a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n + b_n = 0 + a_{11} x_1 + a_{12} x_2 + {} & \ldots & {} + a_{1n} x_n + b_1 = 0 \\ + a_{21} x_1 + a_{22} x_2 + {} & \ldots & {} + a_{2n} x_n + b_2 = 0 \\ + & \vdots & \\ + a_{n1} x_1 + a_{n2} x_2 + {} & \ldots & {} + a_{nn} x_n + b_n = 0 + +.. _linear systems of equations: https://en.wikipedia.org/wiki/System_of_linear_equations +.. _a number of ways to do this: https://docs.sympy.org/latest/guides/solving/index.html These equations can be put into matrix form: @@ -596,16 +776,14 @@ where: a_{21} & a_{22} & \ldots & a_{2n} \\ \ldots & \ldots & \ldots & \ldots \\ a_{n1} & a_{n2} & \ldots & a_{nn} - \end{bmatrix} - + \end{bmatrix}, \bar{x} = \begin{bmatrix} x_1 \\ x_2 \\ \ldots \\ x_n - \end{bmatrix} - + \end{bmatrix}, \bar{b} = \begin{bmatrix} -b_1 \\ @@ -614,15 +792,16 @@ where: -b_n \end{bmatrix} -Finally, :math:`\bar{x}` is found with matrix inversion (if the matrix is +:math:`\bar{x}`, the solution, is found with matrix inversion (if the matrix is invertible): .. math:: \bar{x} = \mathbf{A}^{-1}\bar{b} -Taking the inverse is not computationally efficient, so some form of `Gaussian -elmination`_ should be used to solve the system. +Taking the inverse is not computationally efficient and potentially numerically +inaccurate, so some form of `Gaussian elmination`_ should be used to solve the +system. .. _Gaussian elmination: https://en.wikipedia.org/wiki/Gaussian_elimination @@ -653,16 +832,24 @@ terms that are not linear in :math:`a_1` and :math:`a_2`. .. jupyter-execute:: - b = -exprs.xreplace({a1: 0, a2:0}) + b = -exprs.xreplace({a1: 0, a2: 0}) b -Lastly, the ``LUsolve()`` method performs Gaussian-Elimination to solve the -system: +The ``.inv()`` method can compute the inverse of A to find the solution: + +.. jupyter-execute:: + + A.inv() @ b + +But it is best to use the ``.LUsolve()`` method to perform Gaussian-Elimination +to solve the system, especially as the dimension of :math:`\mathbf{A}` grows: .. jupyter-execute:: A.LUsolve(b) +.. todo:: Add exercise + Simplification ============== @@ -684,7 +871,12 @@ trigonometric simplifications, for example: .. jupyter-execute:: - sm.trigsimp(sm.cos(omega)**2 + sm.sin(omega)**2) + trig_expr = sm.cos(omega)**2 + sm.sin(omega)**2 + trig_expr + +.. jupyter-execute:: + + sm.trigsimp(trig_expr) .. warning:: @@ -695,8 +887,8 @@ trigonometric simplifications, for example: As mentioned earlier, SymPy represents expressions as trees. Symbolic expressions can also be represented as `directed acyclic graphs`_ that contain only one node for each unique expression (unlike SymPy's trees which may have -the same expression in more than one node). These unique expressions, or "common -subexpressions", can be found with the +the same expression in more than one node). These unique expressions, or +"common subexpressions", can be found with the :external:py:func:`~sympy.simplify.cse_main.cse` function. This function will provide a simpler form of the equations that minimizes the number of operations to compute the answer. @@ -766,6 +958,8 @@ else has asked your question), you can do so at the following places: - `SymPy mailing list `_: Ask questions via email. +- `SymPy Github Discussions `_: Ask + questions via Github. - `SymPy Gitter `_: Ask questions in a live chat. - `Stackoverflow From 9fbe92485e8a7f022d0d83ec9695a26900e61e5b Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Tue, 7 Feb 2023 08:16:51 +0100 Subject: [PATCH 2/8] Added link to LU decomposition. --- sympy.rst | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sympy.rst b/sympy.rst index 3ebeee53..fccf8bd4 100644 --- a/sympy.rst +++ b/sympy.rst @@ -841,13 +841,16 @@ The ``.inv()`` method can compute the inverse of A to find the solution: A.inv() @ b -But it is best to use the ``.LUsolve()`` method to perform Gaussian-Elimination -to solve the system, especially as the dimension of :math:`\mathbf{A}` grows: +But it is best to use the ``.LUsolve()`` method to perform `LU decomposition`_ +style Gaussian-Elimination to solve the system, especially as the dimension of +:math:`\mathbf{A}` grows: .. jupyter-execute:: A.LUsolve(b) +.. _LU decomposition: https://en.wikipedia.org/wiki/LU_decomposition + .. todo:: Add exercise Simplification From 0cf78cfbd1e2aefcc75e9c1fbf99c1d5d7f092da Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Tue, 7 Feb 2023 15:49:34 +0100 Subject: [PATCH 3/8] Removed multiequation alignment because pdf wouldnt' build. --- sympy.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sympy.rst b/sympy.rst index fccf8bd4..e3922224 100644 --- a/sympy.rst +++ b/sympy.rst @@ -752,10 +752,10 @@ below. First, you should know you have equations of this form: .. math:: - a_{11} x_1 + a_{12} x_2 + {} & \ldots & {} + a_{1n} x_n + b_1 = 0 \\ - a_{21} x_1 + a_{22} x_2 + {} & \ldots & {} + a_{2n} x_n + b_2 = 0 \\ - & \vdots & \\ - a_{n1} x_1 + a_{n2} x_2 + {} & \ldots & {} + a_{nn} x_n + b_n = 0 + a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n + b_1 = 0 \\ + a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n + b_2 = 0 \\ + \vdots \\ + a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n + b_n = 0 .. _linear systems of equations: https://en.wikipedia.org/wiki/System_of_linear_equations .. _a number of ways to do this: https://docs.sympy.org/latest/guides/solving/index.html From 0073ea0fa4109fecffd4875ade74aaa98ea30cf3 Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Tue, 7 Feb 2023 16:42:08 +0100 Subject: [PATCH 4/8] Completed learning objectives and added two examples. --- sympy.rst | 109 +++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 100 insertions(+), 9 deletions(-) diff --git a/sympy.rst b/sympy.rst index e3922224..d5e20e95 100644 --- a/sympy.rst +++ b/sympy.rst @@ -14,9 +14,12 @@ Learning Objectives After completing this chapter readers will be able to: - Write mathematical expressions with symbols and functions using SymPy. +- Print different forms of expressions and equations with SymPy. - Differentiate mathematical expressions using SymPy. - Evaluate mathematical expressions using SymPy. +- Create matrices and do linear algebra using SymPy. - Solve a linear system of equations with SymPy. +- Simplify mathematical expressions with SymPy. Introduction ============ @@ -740,7 +743,36 @@ respect to a vector (or sequence) of variables. .. _Jacobian: https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant -.. todo:: Add an exercise here. +.. admonition:: Exercise + + Write your own function that produces a Jacobian given a column matrix of + expressions. It should look like:: + + def jacobian(v, x): + """Returns the Jacobian of the vector function v with respect to the + vector of variables x.""" + # fill in your code here + return J_v_x + + Show that it gives the same solution as the above ``.jacobian()`` method. Do + not use the ``.jacobian()`` method in your function. + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + def jacobian(v, x): + """Returns the Jacobian of the vector function v with respect to the + vector of variables x.""" + diffs = [] + for expr in v: + for var in x: + diffs.append(expr.diff(var)) + J_v_x = sm.Matrix(diffs).reshape(len(v), len(x)) + return J_v_x + + jacobian(mat3, mat4) Solving Linear Systems ====================== @@ -841,9 +873,9 @@ The ``.inv()`` method can compute the inverse of A to find the solution: A.inv() @ b -But it is best to use the ``.LUsolve()`` method to perform `LU decomposition`_ -style Gaussian-Elimination to solve the system, especially as the dimension of -:math:`\mathbf{A}` grows: +But it is best to use the ``.LUsolve()`` method to perform an `LU +decomposition`_ Gaussian-Elimination to solve the system, especially as the +dimension of :math:`\mathbf{A}` grows: .. jupyter-execute:: @@ -851,17 +883,74 @@ style Gaussian-Elimination to solve the system, especially as the dimension of .. _LU decomposition: https://en.wikipedia.org/wiki/LU_decomposition -.. todo:: Add exercise +.. admonition:: Exercise + + Solve the following equations for all of the :math:`L`'s and then use + ``lambdify()`` to evaluate the solution for :math:`F_1=13` and + :math:`F_2=32`. + + .. math:: + + -L_1 + L_2 - L_3/\sqrt{2} = & 0 \\ + L_3/\sqrt{2} + L_4 = & F_1 \\ + -L_2 - L_5/\sqrt{2} = & 0 \\ + L_5/\sqrt{2} = & F_2 \\ + L_5/\sqrt{2} + L_6 = & 0 \\ + -L_4 -L_5/\sqrt{2} = & 0 + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + L1, L2, L3, L4, L5, L6, F1, F2 = sm.symbols('L1, L2, L3, L4, L5, L6, F1, F2') + + exprs = sm.Matrix([ + -L1 + L2 - L3/sm.sqrt(2), + L3/sm.sqrt(2) + L4 - F1, + -L2 - L5/sm.sqrt(2), + L5/sm.sqrt(2) - F2, + L5/sm.sqrt(2) + L6, + -L4 -L5/sm.sqrt(2), + ]) + exprs + + .. jupyter-execute:: + + unknowns = sm.Matrix([L1, L2, L3, L4, L5, L6]) + + coef_mat = exprs.jacobian(unknowns) + rhs = exprs.xreplace(dict(zip(unknowns, [0]*6))) + + sol = coef_mat.LUsolve(rhs) + + sm.Eq(unknowns, sol) + + .. jupyter-execute:: + + eval_sol = sm.lambdify((F1, F2), sol) + eval_sol(13, 32) Simplification ============== The above result from :external:py:meth:`~sympy.matrices.matrices.MatrixBase.LUsolve` is a bit -complicated. SymPy has some functionality for automatically simplifying -symbolic expressions. The function -:external:py:func:`~sympy.simplify.simplify.simplify` will attempt to find a -simpler version: +complicated. Reproduced here: + +.. jupyter-execute:: + + a1, a2 = sm.symbols('a1, a2') + exprs = sm.Matrix([ + [a1*sm.sin(f(t))*sm.cos(2*f(t)) + a2 + omega/sm.log(f(t), t) + 100], + [a1*omega**2 + f(t)*a2 + omega + f(t)**3], + ]) + A = exprs.jacobian([a1, a2]) + b = -exprs.xreplace({a1: 0, a2: 0}) + +SymPy has some functionality for automatically simplifying symbolic +expressions. The function :external:py:func:`~sympy.simplify.simplify.simplify` +will attempt to find a simpler version: .. jupyter-execute:: @@ -935,6 +1024,8 @@ intermediate variables. simplified[0] +.. todo:: Add exercise. + Learn more ========== From 51f70a1367073068681148532df5047287dedc9e Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Tue, 7 Feb 2023 22:23:15 +0100 Subject: [PATCH 5/8] Added a lambdify cse exercise to show the speed up. --- sympy.rst | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 54 insertions(+), 4 deletions(-) diff --git a/sympy.rst b/sympy.rst index d5e20e95..ad972521 100644 --- a/sympy.rst +++ b/sympy.rst @@ -947,6 +947,7 @@ complicated. Reproduced here: ]) A = exprs.jacobian([a1, a2]) b = -exprs.xreplace({a1: 0, a2: 0}) + sol = A.LUsolve(b) SymPy has some functionality for automatically simplifying symbolic expressions. The function :external:py:func:`~sympy.simplify.simplify.simplify` @@ -954,7 +955,7 @@ will attempt to find a simpler version: .. jupyter-execute:: - sm.simplify(A.LUsolve(b)) + sm.simplify(sol) But you'll have the best luck at simplifying if you use simplification functions that target the type of expression you have. The @@ -983,13 +984,20 @@ the same expression in more than one node). These unique expressions, or "common subexpressions", can be found with the :external:py:func:`~sympy.simplify.cse_main.cse` function. This function will provide a simpler form of the equations that minimizes the number of operations -to compute the answer. +to compute the answer. We can count the number of basic operations (additions, +multiplies, etc.) using :external:py:func:`~sympy.core.function.count_ops`: + +.. jupyter-execute:: + + sm.count_ops(sol) .. _Directed acyclic graphs: https://en.wikipedia.org/wiki/Directed_acyclic_graph +We can simplify with ``cse()``: + .. jupyter-execute:: - substitutions, simplified = sm.cse(A.LUsolve(b)) + substitutions, simplified = sm.cse(sol) The ``substitutions`` variable contains a list of tuples, where each tuple has a new intermediate variable and the sub-expression it is equal to. @@ -1024,7 +1032,49 @@ intermediate variables. simplified[0] -.. todo:: Add exercise. +We can count the number of operations of the simplified version: + +.. jupyter-execute:: + + num_ops = sm.count_ops(simplified[0]) + for sub in substitutions: + num_ops += sm.count_ops(sub[1]) + num_ops + +.. admonition:: Exercise + + :external:py:func:`~sympy.utilities.lambdify.lambdify` has an optional + argument ``cse=True|False`` that applies common subexpression elimination + internally to simplify the number of operations. Differentiate the + ``base_expr`` with respect to ``x`` 10 times to generate a very long + expression. Create two functions using ``lambdify()``, one with ``cse=True`` + and one with ``cse=False``. Compare how long it takes to numerically + evaluate the resulting functions using the ``%timeit`` magic. + + .. jupyter-execute:: + + a, b, c, x, y, z = sm.symbols('a, b, c, x, y, z') + base_expr = a*sm.sin(x*x + b*sm.cos(x*y) + c*sm.sin(x*z)) + +.. admonition:: Solution + :class: dropdown + + .. jupyter-execute:: + + long_expr = base_expr.diff(x, 10) + + eval_long_expr = sm.lambdify((a, b, c, x, y, z), long_expr) + eval_long_expr_cse = sm.lambdify((a, b, c, x, y, z), long_expr, cse=True) + + .. jupyter-execute:: + + %%timeit + eval_long_expr(1.0, 2.0, 3.0, 4.0, 5.0, 6.0) + + .. jupyter-execute:: + + %%timeit + eval_long_expr_cse(1.0, 2.0, 3.0, 4.0, 5.0, 6.0) Learn more ========== From 32767424ce13b2bd98e1eda70a81f144a484820c Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Wed, 8 Feb 2023 10:56:35 +0100 Subject: [PATCH 6/8] Added warning about LUsolve not ensuring correct solutions and spaced out the solution to last exercise. --- sympy.rst | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/sympy.rst b/sympy.rst index ad972521..2ae56434 100644 --- a/sympy.rst +++ b/sympy.rst @@ -883,6 +883,24 @@ dimension of :math:`\mathbf{A}` grows: .. _LU decomposition: https://en.wikipedia.org/wiki/LU_decomposition +.. warning:: + + This method of solving symbolic linear systems is fast, but it can give + incorrect answers for: + + 1. expressions that are not acutally linear in the variables the Jacobian is + taken with respect to + 2. A matrix entries that would evaluate to zero if simplified or specific + numerical values are provided + + So only use this method if you are sure your equations are linear and if + your A matrix is made up of complex expressions, watch out for ``nan`` + results after lambdifying. :external:py:func:`~sympy.solvers.solvers.solvee` + and :external:py:func:`~sympy.solvers.solveset.linsolve` can also solve + linear systems and they check for linearity and properties of the A matrix. + The cost is that they can be extremely slow for large expressions (which we + will have in this course). + .. admonition:: Exercise Solve the following equations for all of the :math:`L`'s and then use @@ -1059,13 +1077,21 @@ We can count the number of operations of the simplified version: .. admonition:: Solution :class: dropdown + Differentiate 10 times: + .. jupyter-execute:: long_expr = base_expr.diff(x, 10) + Create the numerical functions: + + .. jupyter-execute:: + eval_long_expr = sm.lambdify((a, b, c, x, y, z), long_expr) eval_long_expr_cse = sm.lambdify((a, b, c, x, y, z), long_expr, cse=True) + Now time each function: + .. jupyter-execute:: %%timeit From 77eca9f3da497a11dbe08d06872804ccfed085ba Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Wed, 8 Feb 2023 11:11:50 +0100 Subject: [PATCH 7/8] Spelling and bold A matrix. --- sympy.rst | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/sympy.rst b/sympy.rst index 2ae56434..de162890 100644 --- a/sympy.rst +++ b/sympy.rst @@ -890,16 +890,17 @@ dimension of :math:`\mathbf{A}` grows: 1. expressions that are not acutally linear in the variables the Jacobian is taken with respect to - 2. A matrix entries that would evaluate to zero if simplified or specific - numerical values are provided + 2. :math:`\mathbf{A}` matrix entries that would evaluate to zero if + simplified or specific numerical values are provided So only use this method if you are sure your equations are linear and if - your A matrix is made up of complex expressions, watch out for ``nan`` - results after lambdifying. :external:py:func:`~sympy.solvers.solvers.solvee` - and :external:py:func:`~sympy.solvers.solveset.linsolve` can also solve - linear systems and they check for linearity and properties of the A matrix. - The cost is that they can be extremely slow for large expressions (which we - will have in this course). + your :math:`\mathbf{A}` matrix is made up of complex expressions, watch out + for ``nan`` results after lambdifying. + :external:py:func:`~sympy.solvers.solvers.solve` and + :external:py:func:`~sympy.solvers.solveset.linsolve` can also solve linear + systems and they check for linearity and properties of the A matrix. The + cost is that they can be extremely slow for large expressions (which we will + have in this course). .. admonition:: Exercise From 4a42be4de8590a854b13fb2c9deea95470f2994e Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Wed, 8 Feb 2023 22:06:50 +0100 Subject: [PATCH 8/8] Update sympy.rst --- sympy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympy.rst b/sympy.rst index de162890..0778a148 100644 --- a/sympy.rst +++ b/sympy.rst @@ -40,7 +40,7 @@ and the primary SymPy features needed we will be using. Import and Setup ================ -I will import SymPy as follows through this book: +I will import SymPy as follows throughout this book: .. jupyter-execute::