ex2.ipynb

ex2.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Neutrino Physics, Tutorial Session #2\n",
        "# 2020-05-06"
      ],
      "metadata": {
        "nteract": {
          "transient": {
            "deleting": false
          }
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Today we try to use SymPy to solve our exercise problems."
      ],
      "metadata": {
        "nteract": {
          "transient": {
            "deleting": false
          }
        }
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sympy import *\n",
        "init_printing(use_latex=\"mathjax\")"
      ],
      "outputs": [],
      "execution_count": 1,
      "metadata": {
        "collapsed": true,
        "outputExpanded": false,
        "jupyter": {
          "source_hidden": false,
          "outputs_hidden": false
        },
        "nteract": {
          "transient": {
            "deleting": false
          }
        },
        "execution": {
          "iopub.status.busy": "2020-05-06T13:18:33.789Z",
          "iopub.execute_input": "2020-05-06T13:18:33.792Z",
          "iopub.status.idle": "2020-05-06T13:18:34.649Z",
          "shell.execute_reply": "2020-05-06T13:18:34.651Z"
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "SymPy's `solve` method can be used to solve an equation with respect to a variable. In order to solve systems of equations with SymPy, we need to implement a `focus` method, which is not part of SymPy core, yet.\n",
        "\n",
        "https://stackoverflow.com/a/61547194/2066546\n",
        "https://github.com/sympy/sympy/issues/2720#issuecomment-312437508"
      ],
      "metadata": {
        "nteract": {
          "transient": {
            "deleting": false
          }
        }
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def seek(eqs, do, sol=[], strict=True):\n",
        "    from sympy.solvers.solvers import _invert as f\n",
        "    from sympy.core.compatibility import ordered\n",
        "    while do and eqs:\n",
        "        for x in do:\n",
        "            for e in eqs:\n",
        "                if not isinstance(e, Eq):\n",
        "                    continue\n",
        "                i, d = f(e.lhs - e.rhs, x)\n",
        "                if d != x:\n",
        "                    continue\n",
        "                break\n",
        "            else:\n",
        "                if strict:\n",
        "                    assert None  # no eq could be solved for x\n",
        "                continue\n",
        "            sol.append((d, i))\n",
        "            eqs.remove(e)\n",
        "            break\n",
        "        do.remove(x)\n",
        "        if not strict:\n",
        "            do.extend(i.free_symbols)\n",
        "            do = list(ordered(do))\n",
        "        for _ in range(len(eqs)):\n",
        "            if not isinstance(eqs[_], Eq):\n",
        "                continue\n",
        "            # avoid dividing by zero\n",
        "            ln, ld = eqs[_].lhs.as_numer_denom()\n",
        "            rn, rd = eqs[_].rhs.as_numer_denom()\n",
        "            eqs[_] = Eq(ln*rd, rn*ld).xreplace({x: i})\n",
        "            if eqs[_] == False:\n",
        "                raise ValueError('inconsistency detected')\n",
        "    return sol\n",
        "\n",
        "def focus(eqs, *syms, **kwargs):\n",
        "    \"\"\"Given Equality instances in ``eqs``, solve for symbols in\n",
        "    ``syms`` and resolve as many of the free symbols in the solutions\n",
        "    as possible. When ``evaluate=True`` a dictionary with keys being\n",
        "    ``syms`` is returned, otherwise a list of identified symbols\n",
        "    leading to the desired symbols is given.\n",
        "\n",
        "    Examples\n",
        "    ========\n",
        "    >>> focus((Eq(a, b), Eq(b + 2, c)), a)\n",
        "    {a: c - 2}\n",
        "    >>> focus((Eq(a, b), Eq(b + 2, c)), a, evaluate=False)\n",
        "    [(b, c - 2), (a, b)]\n",
        "    \"\"\"\n",
        "    from sympy.solvers.solvers import _invert as f\n",
        "    from sympy.core.compatibility import ordered\n",
        "    evaluate = kwargs.get('evaluate', True)\n",
        "    assert all(isinstance(i, Eq) for i in eqs)\n",
        "    sol = []\n",
        "    free = Tuple(*eqs).free_symbols\n",
        "    do = set(syms) & free\n",
        "    if not do:\n",
        "        return sol\n",
        "    eqs = list(eqs)\n",
        "    seek(eqs, do, sol)\n",
        "    assert not do\n",
        "    for x, i in sol:\n",
        "        do |= i.free_symbols\n",
        "    do = list(ordered(do))  # make it canonical\n",
        "    seek(eqs, do, sol, strict=False)\n",
        "    if evaluate:\n",
        "        while len(sol) > len(syms):\n",
        "            x, s = sol.pop()\n",
        "            for i in range(len(sol)):\n",
        "                sol[i] = (sol[i][0], sol[i][1].xreplace({x: s}))\n",
        "        for i in reversed(range(1, len(syms))):\n",
        "            x, s = sol[i]\n",
        "            for j in range(i):\n",
        "                y, t = sol[j]\n",
        "                sol[j] = y, f(y - t.xreplace({x: s}), y)[0]\n",
        "    if evaluate:\n",
        "        sol = dict(sol)\n",
        "    else:\n",
        "        sol = list(reversed(sol))\n",
        "    return sol"
      ],
      "outputs": [],
      "execution_count": 2,
      "metadata": {
        "collapsed": true,
        "outputExpanded": false,
        "jupyter": {
          "source_hidden": false,
          "outputs_hidden": false
        },
        "nteract": {
          "transient": {
            "deleting": false
          }
        },
        "execution": {
          "iopub.status.busy": "2020-05-06T13:18:54.939Z",
          "iopub.execute_input": "2020-05-06T13:18:54.942Z",
          "iopub.status.idle": "2020-05-06T13:18:54.945Z",
          "shell.execute_reply": "2020-05-06T13:18:54.948Z"
        }
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": true,
        "outputExpanded": false,
        "jupyter": {
          "source_hidden": false,
          "outputs_hidden": false
        },
        "nteract": {
          "transient": {
            "deleting": false
          }
        }
      }
    }
  ],
  "metadata": {
    "kernel_info": {
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.7.6",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "kernelspec": {
      "argv": [
        "/usr/local/opt/python/bin/python3.7",
        "-m",
        "ipykernel_launcher",
        "-f",
        "{connection_file}"
      ],
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}