{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Riding the wave" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We've reached the third lesson of the module _Riding the wave: Convection problems_, where we explore the numerical solution of conservation laws. We learned in the [first lesson](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb) how to think about conservation laws, starting with the most intuitive case: that of conservation of mass. But there are many physical quantities that can be conserved: energy, for example. Or cars on a road.\n", "\n", "Developing a conservation law for traffic in a one-lane roadway is fun, because we can relate to the insights it gives. We've all experienced traffic slowing down to a crawl when the density of cars gets very high, and stepping on the accelerator pedal with glee when there are no cars on the road!\n", "\n", "In the previous two lessons, we developed a traffic-flow model and explored [different choices of numerical scheme](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_02_convectionSchemes.ipynb). Not everything worked as you expected, and by now you might be realizing how some restrictions come about from the numerical methods themselves, while others still are imposed by the models we use. This lesson will develop a better traffic model, and also show you some impressive SymPy kung fu." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Traffic flow, revisited" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A better flux model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Like you saw in the first lesson, cars obey a typical conservation law,\n", "\n", "\\begin{equation}\n", " \\frac{\\partial \\rho}{\\partial t} + \\frac{\\partial F}{\\partial x} = 0, \\end{equation}\n", "\n", "where $F$ is the flux, $F=\\rho u$—flux equals density times velocity. From our experience on the road, we know that the traffic speed is a function of traffic density, and we proposed as a first approximation a linear relation between the two. Thus,\n", "\n", "\\begin{equation}\n", "F(\\rho) = \\rho \\, u_{max}\\left(1 - \\frac{\\rho}{\\rho_{max}} \\right)\n", "\\end{equation}\n", "\n", "This flux model meets the two requirements, based on a qualitative view of traffic flow, that:\n", "1. $u \\rightarrow u_{max}$ and $F\\rightarrow 0$ when $\\rho \\rightarrow 0$.\n", "2. $u \\rightarrow 0$ as $\\rho \\rightarrow \\rho_{max}$\n", "\n", "However, it leads to some unrealistic or at least improbable results. For example, note that if the traffic speed is a linear function of density, the flux function will be quadratic (see Figure 1). In this case, the maximum flux will occur when $\\rho^{\\star} = \\rho_{max}/2$, corresponding to a traffic speed $u_{max}/2$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![velocity_and_flux](./figures/velocity_and_flux.png)\n", "#### Figure 1. Traffic speed (left) and flux (right) vs. density." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A good question to ask here is: should the maximum flux on a given stretch of road have a strict dependence on the maximum speed that the roadway allows, be it by physical restrictions or posted speed limits? In other words, do we really expect the maximum flux to increase if we allow arbitrarily high speeds? \n", "\n", "Probably not. But there *should* be some ideal traffic speed, $u^{\\star}$, corresponding to an ideal traffic density, $\\rho^{\\star}$, resulting in the maximum traffic flux:\n", "\n", "\\begin{equation}F_{\\rm max} = \\rho^{\\star}u^{\\star}\\end{equation}\n", "\n", "Let us improve the initial flux model by taking this observation into account. One thing that we can try is to introduce a flux model that is cubic in $\\rho$ instead of quadratic:\n", "\n", "\\begin{equation}F(\\rho) = u_{\\rm max}\\rho (1 - A\\rho - B \\rho^2)\\end{equation}\n", "\n", "This new model still meets the first criterion listed above: $F\\rightarrow 0$ when $\\rho \\rightarrow 0$. Moreover, we impose the following conditions:\n", "\n", "* When $\\rho = \\rho_{\\rm max}$ traffic flux goes to zero:\n", "\n", "\\begin{equation}F(\\rho_{\\rm max}) = 0 = u_{\\rm max}\\, \\rho_{\\rm max}(1 - A \\rho_{\\rm max} - B \\rho_{\\rm max}^2)\\end{equation}\n", "\n", "* Based on eq. (3), maximum flux occurs when $\\rho = \\rho^{\\star}$ and $F'(\\rho^{\\star}) = 0$:\n", "\n", "\\begin{equation}\n", "F'(\\rho^{\\star}) = 0 = u_{\\rm max}(1 - 2A\\rho^{\\star} - 3B(\\rho^{\\star})^2)\\end{equation}\n", "\n", "* $u^{\\star}$ is obtained when $\\rho = \\rho^{\\star}$:\n", "\n", "\\begin{equation}\n", "u^{\\star} = u_{\\rm max}(1 - A \\rho^{\\star} - B(\\rho^{\\star})^2)\\end{equation}\n", "\n", "We have three equations and four unknowns $A,B,\\rho^{\\star}, u^{\\star}$. However, in practice, the ideal traffic speed could be obtained for a given road by observations. Similarly to $u_{max}$ and $\\rho_{max}$ it will therefore be taken as a parameter." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving the new flux equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equations $(5)$, $(6)$, and $(7)$ are not incredibly difficult to solve with pen and paper. Instead of following that route, we can use [SymPy](http://sympy.org/en/index.html), the symbolic mathematics library of Python. You used SymPy already in the [Burgers' equation lesson](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_04_1DBurgers.ipynb) of Module 2, _\"Space & Time\"_, and you will learn some new functionalities here.\n", "\n", "We begin by importing SymPy, initializing $\\LaTeX$ printing and defining a set of symbolic variables that we'll use in the calculations. Remember: variables are not defined automatically, you have to tell SymPy that a name will correspond to a symbolic variable by using the keyword `symbols`. This behavior is different from many other symbolic math systems that implicitly construct symbols for you, so you may be perplexed if you've used these other sytems before. The reason for this behavior in SymPy is that the system is fully built on Python, a general-purpose language that needs you to define all objects before using them." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import sympy\n", "sympy.init_printing()\n", "\n", "u_max, u_star, rho_max, rho_star, A, B = sympy.symbols('u_max u_star rho_max rho_star A B')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we used the special character `_` (the under-dash) to create a symbol with a subscript. Here, for example, we've created the symbols $u_{max}$ and $u_{star}$, representing $u_{\\rm max}$ and $u^{\\star}$ ($u$-star) in the equations, and assigned them to the variable names `u_max` and `u_star`. SymPy also allows you to create symbols with a superscript by means of the special character `^`. Be careful, though: SymPy is built on Python, so you denote an exponent in an expression by `**` (this may also be different from other symbolic math systems that you have used in the past).\n", "\n", "Next, use `sympy.Eq()` to define the three equations—corresponding to Equations $(5)$, $(6)$ and $(7)$, above—in terms of symbolic variables. The function `sympy.Eq()` creates an equality between two SymPy objects, passed as arguments separated by a comma. We need to remember that the equal sign in Python is used for variable assignment, and for that reason it cannot be used to build a symbolic equation with SymPy. That is why we need `sympy.Eq()` to create symbolic equalities. But the equal sign _is_ used to _assign_ an equation to a name: here, `eq1`, `eq2`, `eq3` are names for the symbolic equalities we create with `sympy.Eq()`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "eq1 = sympy.Eq( 0, u_max*rho_max*(1 - A*rho_max-B*rho_max**2) )\n", "eq2 = sympy.Eq( 0, u_max*(1 - 2*A*rho_star-3*B*rho_star**2) )\n", "eq3 = sympy.Eq( u_star, u_max*(1 - A*rho_star - B*rho_star**2) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check this out: you can display these equations in pretty typeset mathematics just by executing their name in a code cell:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAAaCAMAAAAuYg+iAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAABGhJREFUaAXtWuHW5BAMpVqqSr3/024S\nlOm0hn79zs7Zs350FLniiggdxv6nAwNKCnMo+gde+TeMCYhVs/0hm+KH8o+LK/k4ZD+g9Rtjcu4X\nfJH4iqEUGo1L8fLXssYBt5v/af9a/RThUXnNa3Dm4WVmVUxnnerprLSnzETL58Mgxh7BX2lr1irs\n5KrV3ZVC+lVCcuv7lNrlNh2jjppMBDuBm1PzbbTuUV0I6KTVab3yD3PLjA8dTv64fO10kwwj9Zy8\nCUfDpQerj+x0uA8XruDnrpNc62Z9LXlVw30wWOkPW6gdgJQrqUO5ePFUozIycTt6iDUW8i3C35yr\nQ2e3X62vBWBidEnp2z0cBHXscIgcp2qrwREDvU1JHuaF7dyyGWgP06beVsYJtlGNE2D640NR485w\nVnLbrMfJGPaiOew247qHJwRrVg/pAW7dxFTi9jgDuxJ7ZpDKNZG2cY1oW88yroaUgDbkyW/WY9f8\nJBPdrV1cspc7sBW71e6dW43zllK5vzgw8ia6RsEG3Hl0TwTuKnuVgt1G7tw263HCaC7ifoIoQfvd\nYd6CrXAL3qHZbiWunSbXAUfXGSeFHnk09VyNW1yhMrnFdj2qHUZ3q9Lyvwf7DLcjxYF8t56K4nAg\nx2YmsVFpmqsit8YVKe5uUpOFBafVoQeBnwFiRXS3bApbWi8sk6TnPNNP9qd5L4OcSf72Q9jBSYe0\nuWZOTnPkDYLEaf1JYbLbdFzC39DMSMyLGCp16YEAKJxS7jdFt0wHa+mGNYSpMahQKoc4L9yylXwq\nx3CM0oanlZRyzBmc54xeEaIF3FIhmgsSYyQBf0MkQd4AJLDB59aEkrgNkOUz7NkmHh669ChhXvMp\numVzsNsC9rVh/a3uExjsZ5DyNn3BLdFFa11sy8a3xQq1ELl8UIRhxebkRhEN2RhIwFppaE3qZwUO\no0nmG7nt0uMAVbymBSiiv82wRaPP2Qq3yIlZwWLHOVvoKWJYQws4DsMN2jg6U4cGuwGGRNsHImdj\nyMpWyFsgWDa1pg6v4luTrk08zVqfHqdDocLobkU8+GbY5oVGMBVuIb4FLMe3aQ9FLtThXkCr4JM5\nMEn3PDPa7QyF8QorXf4w7hTfBk33959bU5fn5zLjfLDXYYEMKNmnx8VgjFv8ipvQOkV3VsC2LjTC\nPnALiH5xgSU8l0FSwUtS/uKBzjPqQUwiY3ZBh4r+CuwV5YhGAiBfbKkQef/QmkTq9wnUBB6deiSx\nT78ZtmVZZrQDt7mC0X1C8V7Jkj+K9WiuyJjmcEpg6wi3aLj6VwanMSx5Ta2t226L7uvxqtXhrYRt\nXGiEYFMEcMADc2o+OJWxKn0fcMAvnwY4L3KNe5hl4yS3Sb5R29z6w/1tUP6+Hm+DLwtK2JZlWcpe\n5KP/vKgti7cQqZRFj+fr3x1Cd7+kxwts60KrErBvPdVWofLS9htkG5s0fS/7JT0K2OaFVh3Wl30v\nY9/2cbRKXr3y+4byFf9PqJPWWEvbzh8a9S1/S2d49wAAAABJRU5ErkJggg==\n", "text/latex": [ "$$0 = \\rho_{max} u_{max} \\left(- A \\rho_{max} - B \\rho_{max}^{2} + 1\\right)$$" ], "text/plain": [ " ⎛ 2 ⎞\n", "0 = ρ_max⋅u_max⋅⎝-A⋅ρ_max - B⋅ρ_max + 1⎠" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq1" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAAaCAMAAAD2bquoAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAABPxJREFUaAXlWdu2oyAMBUUQEeT/v3aS\nQBAvtVh7zppZw0NVhJ1kE5JghfhPm9G9+wdMl3+pkkCeGf0jAvtHs9smG9027rdH+bgIocdHYn/B\nNjU90vDnJrsA/C3xmQBrns1/P9vK92NwhPupveBNbmd62OGst73PZf+VXder9mntI918OlbpELYC\nh3A68LNO1a/wvY6zhhbm40r66UOjlc2KDQQ6QIwy44dYlyZaFrQZpbQXLsxV8Dbxi/ypwQs1MDcu\nJh2GuN9sfvjIZKftyPteogPSjzg3dWP3/YcZgsyxWUzJKlbBW8/njnqc29DTRwgGS6ZNyJgcT8dd\nJvMdmN4AJ0S/CS7KOM38qQhOMFEU6ONHq3GpgI+nxUuk7oCyU+tVYI2468HV4KoVj7ZZhy7zyMDe\nQmAEChua3jEvCn9iBGrTwpiDf++QXTF39+L1Y3/OykCxIhSBToqaP2e+sJISnRDbmNxczaUSIHw3\nR2iP+QsDrBTRmy9J6MnvIi2OW+5stMv6aixIgNsVMkWnTbi/VHuNQ468Ofz5KfCi3Me/8D8bjvxZ\nXBlunABULzqM8baKWnudD8/hIisYdhBhILDrwl/A4FWoPUA2dRg9sWQZB8i+tggTH+Bf8Ac7uc3/\n4Bg2Ipf002QEDLrib+KVEbiJNEcnjdvsbSx5q4CyU/LhHP4Mb9VP8L/BH5wW0UUc2/nWAhyQ+XOh\najmjrOleW3KQFKMVlWmyeCOicCWC97t2hpyGQHqnDZvDnxhSGjnD32HWj1CnQhtHuiQN8fWaP+DO\ncfx7k8xp58rzjFrLrO7Z//gMgNf0WpaSwGns7XN1kfA5Z9JYLvIr3PUW53Jbe/Eu0Epz9SdsWpIT\n/O207ZMjbIup2pi1ltjwJ2baSbKUEwvW69zWAo52LpCIMJAhMYtBKZTkqUwLXqvsyfxtlYIniWvl\naHJKgS4X0Cm8jhS9lIEtqIwOhMgi4IGlHmCpw6VMG2hFuPoTY1r4Gp+nX/h3GnK9fwXkEGhrrnzB\nH+kDJIIb98u0yGXyvZmIA9kZAvH9EvRSCoUKk3VNV0+uTic4dsjMHy1SCn9Lb6SGZQoa+WMRILrL\nUreY5clHUmAkh2NP7nP8q/B5wqV/06AL/tBshycpNa6exsjb6wzEeSARbJIOnRU3SEDHWwBEoxPD\niNG5tah6Uf/5QN6NZjo+wSej03abiF6IX1AYImIlAkTPSliuReDNsQ04waTzRw5/fT6N1Pjg2o3+\nfcEf1H8gK8hlKAHpqFDqkcHIpbP00VYCW7RuZNwIZuYPLdvFfHH+oNI1gse5gL/QugluQAOoekEV\nog9KwgkjdkJcRaDo66a6QXe42i5MccbQPw8cdiv8dv/e8Qd4cQpJSTx/QDNV0HqpHEU7TxsW2UJD\n/IThEEML+B1OJGJXhPPz7/p+f4fhNZvqlIMIQYhKVSJavoipV/ZU+KLdv3f8VUrT+bd6br1Ft0ND\nrIS6WsCOMiPu61nAGQV7uK1VCvdcXyk8pSEj8NglMRAoWUTeztcor99W+Df828Mynrf06eX83UUv\nfVUOYJwcOoztFvMGfjTSy6Ar+pDRO62uLjGLUXAdALmIEI8+aNf4T/y72JRDTXn++k3LbluFLreq\ny3Ve612N/w3/3kb7Vi3ujLvpLi83yh2ZF2Mr/G/498///yF+4T+qC75+9tWv2Pa3/v/7BWop0v8B\nbIIsBuAt7IQAAAAASUVORK5CYII=\n", "text/latex": [ "$$0 = u_{max} \\left(- 2 A \\rho_{star} - 3 B \\rho_{star}^{2} + 1\\right)$$" ], "text/plain": [ " ⎛ 2 ⎞\n", "0 = u_max⋅⎝-2⋅A⋅ρ_star - 3⋅B⋅ρ_star + 1⎠" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq2" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAAaCAMAAADrEJUzAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAHi5OXp9+7w5u\nj7/fPs+vjb133AAAAAlwSFlzAAAOxAAADsQBlSsOGwAABFJJREFUaAXtWduypCAMBAVEEeX/v3aT\nQEYcFRXdOqe2lodRGdKGJjdQiP9NeNPLf4IG9fPTACb92Dxjs38m/o60N+/gPEBpQiuEGR8ggOgv\nmIew07M5vCEtHZDZhodQ2j8EeC6uVRFDvuo9jU9t55163um80yUfmvadd+2Pld1+P/fOju/euPYm\ndAaa6zZL2Ey28g1WJ8F5A1qJWCumWZN9AB9eJVPIEN83hy+fbOY6LqXRI8cH9dOm2UG0KjTTnVhu\nQXbvLxWi9ZiwTnvNIMQ1u+rXccd6aZhMGx4WBHsa3+hrQrEu6q1jTW+AFobq9L4hkZqGNhpiKfB5\noZn1KoDEh0wxItHyJUbv4/RFrqQSOZnS17liRtEYPdF2XEMQpuwCtOdkOshhrdJId/vUoypwysUd\naDUswW0w3j1c9RQym8mlZbmPWbJM7YTtxYBhXl+PnxpXkhunkAoc4Ur5xUOSMB8yHfjQ0/VWYYZc\nrgPHvQrMEpng8LCdG5EQ+slc4u5tDU6RTPQ7w8HNoGP6D7V3lYvjU8j0yaVrME/IhH0p6ihZ7aTn\n7brrAKc4bSZTuqylnGQ0WVEM+JYqQ7Ui81jDPTjUI4VMMVMe2sMsaGtIxXGkS5aGlgREd+Tgap1Y\na+r5PZyCdvAXkyl4a4LXKCIN3vephonacTKmESUNUZTbogFXmULTouxgLmO3d5IANSZ+77Mi5ItM\ncnBggkZYD1HeeuMod9o0NXjIAFrcRnBbKsUcBwQwU0IdFtViILwSMnV/yNzqHpOrTFV7DOgjhdgr\nGm7xoIerTDGS4eSYPP7Y2tOIEzcXglYfmKCc3ntlgAdncMpq8JCihOjbafDTZzkOyMxwQKBV7dT0\nSYqBmr51puXCpHBUwwaayKR1iiGzvaIhk7O6smH3MWZmmDysZO1xTIlMoqoDFhtgAskcwaJwLzAi\ncy0QaSA5SSU7K/RZmbfggABuBzAMO7RtBsJXjFJ+SrrDOlPysUMg4qODTrRJqdcwhcw+7ilzTPCg\nc3+EeQAdSNKqLW6OdaZQzqt20HTYPIQJh8clGkH7eKyjuABa4Xw95DgoQCC0KgtQQk6SBzsg6UK0\nyGGCGyhkVOhVmw4SKjWUbgodJo9ujqErx7xu7d9kAl6YHK1y3AHFyNagLUorwTlBey2sFRhawJCg\n/+SkLHFDEZJwSAAJbSYMxQsQIsfgTDIne/OEizWwTNH7mYYMCNcM84Y/fpOZAW725iOoPMDsgbze\nCvBtP0bvR0bvNDRIXAGtoJj/AMGJG+y2sCe1k1MjHkbhLT68p2HmbNetvQFLO2ibUyNMLRTkZkP2\nifkCo8knZRwAfXeTgANC1TxgItMRyM6mnc3CpTg7z4y4eQ38loY55jvW/uPnmWS93yuxeW7XNfDm\n/5qOHPMVaz+vBGrUvCVzzeaPfevWy1aDM8xXrP0XfAP6FV/1ViRXPvyGr5MQT+/mtsrZ/mUxSgR/\nAMgmKAtCgcljAAAAAElFTkSuQmCC\n", "text/latex": [ "$$u_{star} = u_{max} \\left(- A \\rho_{star} - B \\rho_{star}^{2} + 1\\right)$$" ], "text/plain": [ " ⎛ 2 ⎞\n", "u_star = u_max⋅⎝-A⋅ρ_star - B⋅ρ_star + 1⎠" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As stressed above, we have three equations with three unknowns (assuming we have some observed value of $u_{star}$, corresponding to the ideal traffic speed)—there must be a solution!\n", "\n", "To eliminate the term with $B$ in `eq2`, leaving it only in terms of $A$ and $\\rho^{\\star}$, we might subtract `3*eq3`. But this will not work in SymPy if you attempt equation subtraction using the equation names. Just like we couldn't use the equal sign (which in Python means assignment) to create a symbolic equality (needing `sympy.Eq` instead), we can't use mathematical operators directly on the SymPy equations, because they are really _equalities_. What does it mean to subtract two _equalities_? It doesn't make sense. Try it ..." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAAaCAMAAACTtSQmAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAAB4BJREFUeAHlW23DsyoItjKtLOv//9oD\nvptWrtzu7Tl+KHMKeEGI2Aj5nxbOOvE/nfqfTftHMG++1DDAYPmwPNJe92j0hwd/gx5+BHPOPqyb\nQnbLthLChsLe+W7fOrectN8g649gTqccgGmbeJfPWrgpO55Cgs2u2671xceRvzjgz7pf6KEu/B/B\nvGnbjr4Fz7HJkaVMyphhL3P9brbRzpPv2DYzKHLOSDL29zjQUY8TBX46FOYetxqj8npwlKvCTz6B\neQ/rNx/eYbRidrAEFcoWIuQcBJN8q2iztF8I7SczH7FpA+u3xCkutk8g23VVsHGw/rnPvAgxhViY\n+LfPPeX14PhXhR+ovh/zRnmL0fgON5EalTzREbdldAucFJuzxn1PhG6DQGM1pkqaTdsV2/aR9dIX\nvqZdFLhQLpi1WQ3dmZixMGc93/lbXg+OY1X4gepzzGPISYL5pFbIbivUoJtpQWWGoDEt24ZGKzfn\naDsqrRWkvV9u4cjVuY5RMSOkNbbryC0tgOuezipsb+3OZqmfwwGBWJiDTm9vzuvBsq0LP1B9jnkC\nOYkx1w6Ip2unnZO+C2dicfvJ02LsZdelV3GIdAxFAwbs+whe4e1p0NliGbQ7p7PbDWr6ywibMzDb\ngpIA6PAjQ+SCj2hpYW4AmCX4Op0DPRjq9eF/jnkCubdZxJxbm917kxiwtRmxw/rKIt4FlhhTg6fB\nUQK6rTNg0jIuX3899uSl3SbpGGGZpH0RNH0xb1Ae26ws2sUpYe4AuJ8VPt+gc6oHUh1+E84+wfzM\nZgHzxGZH1KYtNsqlHWlxnzQGUWgO0ajtLP/JrSMkHDZHzNmsxGDUmXNErviBs8nu6Zqth6zB6JiR\nG/QTAL2fHS2fY9mMMDcArKaIMz28Af4KmCeQB34WME9sNgs/HKMMaL/qku2RaZQnGp3s20DQ2TEb\nbTJcwi/jlAyvuImOk/bVJrTi1qXeoZ8A6G3W12L+0ZMS5haAERn9cIfOmR7eAH8FzBPIA5sFzMts\nFk460RUKa1sZNNMmg5WQQTFHuX4ny0blCHVkQlUatXFeN6UZteQo6w6QllDBgAmtSK9D6xfpQx4Z\nyjCom5YQyXtL9TXARvUyl+jEGoW5BWA0Wf1wh4612ZyED+DPkUMhH2F+AHmMubDx7MUmWkUFTX5X\nlcEWmyxW9iwK77pr4/YugmFrZzJRmr7dd+q+Z1lUHGuL7m2vUr1dNlMIW1nFOUffjsjchaI94naN\nc2+G3lJ9DUarXuYSE1PCZAE8m1xMwz1l6bhfcxWrh4yEj+DPT/gZ5geQxzZLZrVKNy5ts+K5kS0+\nWaWiAsDLqy4HT9TmsIpa4aHB90MoSnobJMyhgtIHGcKYouSwKaQvdIZAqrfAZgrJoF+2DP1wbL6e\nLFTeUn0tPzQUJgfg6eSqKeJQD7D5VXJXhN9nZ59gnkC+s1m9j/Bx+gFUygYAd1wiIRuFGSNI9WpV\nUeM+8R5kqo6wWpRLV6e31vEam1V6NeEs5RCSUs6komhZwIPlmjeTZVNZrUE5VuuxOxPPhvTt8EtX\nlwDoLdXXLLX4HgoTAgjzup5cNUUc6QGXBi1vPfh9dvYJ5gnkO5sVeIpKB+9RY9jt0wzGugDucOvW\naW3Waen4pIy2abky/KVbJVtdIvTou6lFKi+O/YTNFWlD08vKpEx67XjD4NWQDG3WsgDWreFqBdvf\ne5SJ6xyXCa06cyoW0rfDTl2d6pQA6C31Mm8QCBMASG5PDgQK6BQqwvsjO2lzfwf8Npx9hHkCeWCz\nCnMum7V3AeZuWu6xkbxZ2xG+gRSNwEACA0aJ7+kKSznDAAPMeRDCJz0P8oIqNbrBqy0kXqG0E1RA\nAki+gyjKZAnsoSDbjRThYlkA65mS0eZacWxSaNuzFhUr5LTNuCeae+NOQvpFrk4R3wEI9LZJaiEv\n87NOGJicAxAmdXdyIJCnU6qIAz28A/5KmO8gh51RgjkPFvTEBmyDigQW7VjBQpV/UgY1gP7Mh0Ox\n0zo/f7F0gzuGy8a84JhhArvTbIhn0bgcWTBuV6VH8wnpl7u6PYABt4JzMC9MAOCDyWG8BEp4RRGl\negjheSRhgBBm9J1OizE/gbzw7DEUQdXRQtF4lgm3ZLi9Af+Kv2Cjqugh5+fcuk94VeGmbhBUQPSh\nKFIasLj4qi6kltYD+i+4ugXEyJfr7w3y46pNrlARhXoI4KkmYZTRL15ejiGHj6sO1ZEH27Sie0W4\nxgaOdwis1nyAwAC+eIPzXWyxxWdhbcvpPcz+DuBvW80GNmuWhQkVTqkc/xjSf+LHHYfr77pc16hS\nbXKFiijTQwhPNQlDon+JufrqXYJBNX2L+6MR91740SpbexaYLFrxK2UNsr+4aVbBcg+UHQty8cH9\nObeQfhVHcv39bF6gWpMrVUSZHkJ4akkIH4X6dfdPMc9rItP64kp+z+tn+B40BfRrOJI4ej/g+RXN\nZXoI4KkndUD0NzB/5Bbr4ZahVMOR/DP/B8vg846mH8H8G/7v+Q74keYvze2XZD3T10fm8Q3/qz8D\n4cFvYeT+gMxnhv4jelCY/webUlnYmpFtDQAAAABJRU5ErkJggg==\n", "text/latex": [ "$$0 = u_{max} \\left(- 2 A \\rho_{star} - 3 B \\rho_{star}^{2} + 1\\right) - 3 u_{star} = u_{max} \\left(- A \\rho_{star} - B \\rho_{star}^{2} + 1\\right)$$" ], "text/plain": [ " ⎛ 2 ⎞ ⎛ \n", "0 = u_max⋅⎝-2⋅A⋅ρ_star - 3⋅B⋅ρ_star + 1⎠ - 3⋅u_star = u_max⋅⎝-A⋅ρ_star - B⋅ρ_\n", "\n", " 2 ⎞\n", "star + 1⎠" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq2 - 3*eq3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See? SymPy just printed out what you suggested but it did not manipulate algebraically the left and right sides of `eq2` like we were aiming for. What we _can_ do is create a _new_ equation, perform the left-hand side (LHS) and right-hand side (RHS) operations separately and then recombine them into our desired result. For this, it is helpful to know that there are built-in properties of a SymPy equality for the left-hand side and the right-hand side of the equality. Check it out:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqYAAAAaCAMAAABSMm6KAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAHr8+Dm6vLs9O\nn+9eft+PMnn1HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAABvpJREFUeAHtW9m2pCoMZXIW9f+/thMm\no4WIitW1+l4ejoqwSTYxJFCHsf/LVQakrvjVPv/V9v8+V+OP2gLYqJyaZ3ZXPeue1/sHCCzA1Veo\nyiM01krqWO3fr2uWjjE9PRPkG8p9Y4wTFkpw9QNqJLRUbeIlecXf+toa6QoZDG+5ADPtll3t1cde\nXu1xtX2awLK0vcpVnCpVCVGpq6QUb9+PMUil98LVItbuZh3VvdLLoKGIISJJX98dobcd+bk7psLc\nGC1OoAcqSht7hSuVokrVDVN1+7ftlA+eUHpVugFvNpDAUC4FzXSrO18sT/Xy4fmae/xw3U/eDdcR\n46eqwiQ8mog4gX6EorQBaHGuzqiqFlhFOzdDXqvvX3v3KW1H7jGrUgvxRHqI2vO2W+7TVvdxsZak\nl32U3NSZX3G1jUiU5Nqb6UiUiMq3FSbaJFUZJ9D3KEobgBbn6owqOUDoVfpj8+zkX1GKz7IsaKdi\nCe60UsLP+2fryzVb3XszGGOzM9cA18wwL+EpdaP3Bs6CmapViTjCVph4m0RtnEDXoSxtAPoGV+dU\njehSGQ/WkKAj49UNnMaZyA69ntCNibAK8xEe1jZcZnq5tcvnndUd6ifr79QQkjmL3/SQW4GlZpSE\nmbJp62njaLcn4oBAO0p52t7gKpjpIVUCl9xu7NEZdE9X1Ts4FTG+jwmcgkQg3xxsls1aiudfltEd\nxnThVtMKb/sWnw8LlOdmKnKysNsTkSSwOG2vcLWaaZQqqVtMS1TFZrz2Z0EUMaMeZ9AXF13ewUnu\nS0rj6nFUCbmNDmYqMKgOFkykunDrdMce41JDnt+HwdgN/JQ37U9zvwcTkSawOG2vcLWa6QFVqm8b\nBmcYE1qa+XNhpvdNb+GIxBy2IbtCl6Z95KhxbZbBaPdi5D4b3bGxC7ekd5x38FNmus7CsWj3JyJF\nYHnaXuFqJQjuuCDFH09CMq3gSBDnnHszOCYz/eYWjmM5JtyawureuDubpiizvTnmmmkM2eqBuuOd\nC7dYbcPki/iwvwtlmszFSmjgN9zbAYFj08r98XNgXt6eiASBL9D2CldbqvxhC149b5CYQCprlvtx\nl8zc2DGM4oShojfeGXwKN4bMg2t8W7ntIiunzzgtaErWT2QviNHdh6boVQ0tMXzfI3LlBr/HbEtK\nYnlb7n1H08r98XX2ejwR23b7p2MCX6DN75qW5SpOlVWU27xW4OSb5R5sjLDM0qcnHR7b+BK2lGI4\ne1p3z57lXTVEjLgTxI1ENovhbn/fhtATDRbSsn5AM6p72Alkk/1OI/ifCPuaB4s+FaYkgaw4beuu\naVmuUmbaLGb7ZUIXYvwUUGTWLCUhiVZSC7Ppo5zjhQfqK+JmSnGgA27rwO6tnVIPhFeDbKuPzLQx\n+5XmLNf7fmemZipdaJol696kGKO6+9CUVS42pfi+a8pbmzYPzJQKU5BAk3aibOVoC2F8Ya5SZspq\ntB9pTqEGsM8GKEIz7So5arAwodGYxlma7Kvq2lm21NsiA/tCcKBDN3ZtU7leHqipOqG7sEF5lKg2\nwrhqbMf9fo79ruy60xorvi0r0d2HppU7j6P4Xr9zb50y04P01YMfTMRDAl+h7SWuVjONUKXmWs9o\nWGCMQo7d3JvfAENaAZvCEAigTXawtGrIt/nIB8V6v6+IfWJlxYEOePiCga9Af+yB8EuYOLfrEUIc\nbPuZLcsFHAEX+BfK3MINhKuwDz527pz8tqxBdy7aZcCUZqjdskHxYS3IWFlQur2ZAuDSCvMtsehm\nIHZyJQhDJ+IhgW/Q9hJXp1Qp6VdfszQ3xlnOS4umax3IBETbH+CMYXPIsxu7EhzsYECMva9Au6A3\neYgSGwJDaGdRsON/W9ZV990gFD/fW+/NlIKen0KtwrxGIFXrAW1UL5N6h7kowtU5VU4Arjgs1fBd\n90wphkkFOD94l/612FZ684Qd0FSbFjOzFchUmg/CdkoeSUdwTeho64vJSoch+CzfWzdA2UE5PdM/\n6FeWQKJWOdoIaBGu8qmawFPNYEVAEuQvsNLLyYSsxiceEBqtxg7IdD/C6VQAgl+dwcEs1viy7o76\nmuSV7u4Wk5WMSPGfeOsV8vQXUmvTzV1RAqlaxWijoEW4yqcKU3gTTNbwi0/wfJjxYHy2Jj0bLg8f\nTAcBpjrWM6ZiDkjVuqs1sVLYHDrEiL3oyO5uKVnpOBS/jNs5/b0pHX69L0sgVasYbRS0CFc3qVpJ\ne/EOXe6Fcry6XgBJNCX4RdzO+S5BQpisV1kEErWyQLMaEdASXL1PVZZW8UaXHXUc5oXaIm4n/g8+\nJaX9CQJLcPU+VU9o/+3/KHyiGfT9hnLfGOMhDTndf12NH/g38xwab7WhYfgtgJxO/waBhqo/1KhR\ndFso74oAAAAASUVORK5CYII=\n", "text/latex": [ "$$- 3 u_{star} = u_{max} \\left(- 2 A \\rho_{star} - 3 B \\rho_{star}^{2} + 1\\right) - 3 u_{max} \\left(- A \\rho_{star} - B \\rho_{star}^{2} + 1\\right)$$" ], "text/plain": [ " ⎛ 2 ⎞ ⎛ \n", "-3⋅u_star = u_max⋅⎝-2⋅A⋅ρ_star - 3⋅B⋅ρ_star + 1⎠ - 3⋅u_max⋅⎝-A⋅ρ_star - B⋅ρ_s\n", "\n", " 2 ⎞\n", "tar + 1⎠" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq4 = sympy.Eq(eq2.lhs - 3*eq3.lhs, eq2.rhs - 3*eq3.rhs)\n", "eq4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That still needs a little work. SymPy offers several methods to help reduce expressions. You can use [`simplify()`](http://docs.sympy.org/latest/modules/simplify/simplify.html?highlight=simplify#sympy.simplify.simplify.simplify) to attempt to make the expression \"simpler,\" but you can imagine that the quality of an expression being simple is not well defined. One may expect the simpler expression to be shorter, maybe; but not always. The SymPy `simplify()` function applies several strategies, heuristically, to give you an equivalent reduced expression. But some expressions are uncooperative and `simplify()` gives up, returning the argument unchanged. Let's see what it can do with our expression `eq4`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAAVCAMAAABRy6MiAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAHr8+Dm6vLs9O\nn+9eft+PMnn1HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAAA4FJREFUWAnNV9mCnCAQ5BIUQfj/r01X\nI4gJM3FHN1keRuUoqvqCEeKHtuW7ebnv3uAOfurZyekO1Iu1IQwHzKTUZIZD/64znpxv1XM7m1jl\nJT1ANdYJY+f/q9+d9Op8+hyQvt5lohNSrYgtMxI5ZYq0JfvriN8w056cH9f1sT28hPC8ATANckqv\ntPWT9v6AuTzJnYzKH4CMl+QM/SrD/Y6NMJgXEAJC9uVnMOtq15dxJttBy0Bsj2+pb+Wl3bBcZc57\nzoEDur0pxP4SfMTj5Io25frL13Es9q2N3lNhi54Utbrvla1oUoPoFzrOKDVmEglP/ypEwObcfO5a\nrR4f4MzdmaTJT7HJB+Hb/qDcLrp9OvPfv4yfnQhSbJDAP8Np1zo/wGlqaQcwjHm3R5zpS/fD1zj8\nNmveXRMVnQJdQ1lAo9JoZNlH1q3LiBidFvvQ+PECZzy59Hb6oo8x+lyywawwQ+iGxTs+Q2kE4Gtg\nknyhu9Y4KYjmsA9cKtuAvJ4Kbc0Ipw2OXg59MoLdtMsvXHzP6C2fTll3xQkt4yH/3OSK+KLSSPbm\nsCfuHBNGU8UxOiquvGaHow9dQ4aykjzV2nF09zi0ALXb1FUVCE9Gxu5rI1uSU+73nlKGNiZ9hQ+w\n/mh82ZVMOvYnDE90meVvcADbnLgj8pZJh0jMVQTJkLQHiWmZk56b/hfyOxxasIRldtO+qgK5aVFx\n4a0JVpVYR2QyKbHLZzuW1L/Epyw+/zouI+VW3/Y5pliI0XzrW4mFI+4gQ8clHcD0xPBC0iMlkAxy\nNcL/7SA+cGgBLhzILAVlFQhbbFLWQpz2F1m9U3wimVU5Fj7l4xSHZzH0fEQokeFmko0JhMnHSocl\n+QmKU57RV1JtQ2FgK4ZaRTD/VetxsIBB2I4H0I5cIEoYSJWL19NML5SwdBkLi+V9P+azlqO55E8e\nOc7omoScoQ7qpZEUssTAC2MEag85i/qLEV7prv0HDi+ACdyMknIAAbkUGSxi29TV7YkytGfDPT47\nIra81jbaNhFfkktpQxGvt5ITsMFXGoTBZj7QdagB0SWfbpjoKS1y9NWv+uTULx+P8OH8q+Bvnyhq\nnKyW/jCSp1CpcBLUYvV2cTfICxSZINiEEroDGRsXG5t6qgyDsOxvIE/w0bXUdPx+yKsbxOXSn/f3\neRo7MPF92GcQ3CD8KdwebBHqfwFBjR2JVmSougAAAABJRU5ErkJggg==\n", "text/latex": [ "$$- 3 u_{star} = u_{max} \\left(A \\rho_{star} - 2\\right)$$" ], "text/plain": [ "-3⋅u_star = u_max⋅(A⋅ρ_star - 2)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq4.simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That is actually a useful result. We see that `eq4` allows solving for $\\rho^{\\star}$ in terms of $A$, for example (remember that $u_{max}$ is known: the road's maximum velocity). Then we could try to manipulate `eq1` to solve for $B$ in terms of $A$, as well, so that we can substitute back in `eq2` leaving everything in terms of $A$. We can do all that without the help of `simplify()`, but using it interactively to reason about the long expression helped us to decide on a course of action.\n", "\n", "Another SymPy function that can help us examine complicated expressions is [`expand()`](http://docs.sympy.org/latest/modules/core.html?highlight=expand#sympy.core.function.expand). Its purpose is to expand bracketed factors in expressions and group powers of symbols. Let's see what that does here. Notice first that `eq4` hasn't changed; we just printed the result of applying `simplify()` to it." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqYAAAAaCAMAAABSMm6KAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAHr8+Dm6vLs9O\nn+9eft+PMnn1HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAABvpJREFUeAHtW9m2pCoMZXIW9f+/thMm\no4WIitW1+l4ejoqwSTYxJFCHsf/LVQakrvjVPv/V9v8+V+OP2gLYqJyaZ3ZXPeue1/sHCCzA1Veo\nyiM01krqWO3fr2uWjjE9PRPkG8p9Y4wTFkpw9QNqJLRUbeIlecXf+toa6QoZDG+5ADPtll3t1cde\nXu1xtX2awLK0vcpVnCpVCVGpq6QUb9+PMUil98LVItbuZh3VvdLLoKGIISJJX98dobcd+bk7psLc\nGC1OoAcqSht7hSuVokrVDVN1+7ftlA+eUHpVugFvNpDAUC4FzXSrO18sT/Xy4fmae/xw3U/eDdcR\n46eqwiQ8mog4gX6EorQBaHGuzqiqFlhFOzdDXqvvX3v3KW1H7jGrUgvxRHqI2vO2W+7TVvdxsZak\nl32U3NSZX3G1jUiU5Nqb6UiUiMq3FSbaJFUZJ9D3KEobgBbn6owqOUDoVfpj8+zkX1GKz7IsaKdi\nCe60UsLP+2fryzVb3XszGGOzM9cA18wwL+EpdaP3Bs6CmapViTjCVph4m0RtnEDXoSxtAPoGV+dU\njehSGQ/WkKAj49UNnMaZyA69ntCNibAK8xEe1jZcZnq5tcvnndUd6ifr79QQkjmL3/SQW4GlZpSE\nmbJp62njaLcn4oBAO0p52t7gKpjpIVUCl9xu7NEZdE9X1Ts4FTG+jwmcgkQg3xxsls1aiudfltEd\nxnThVtMKb/sWnw8LlOdmKnKysNsTkSSwOG2vcLWaaZQqqVtMS1TFZrz2Z0EUMaMeZ9AXF13ewUnu\nS0rj6nFUCbmNDmYqMKgOFkykunDrdMce41JDnt+HwdgN/JQ37U9zvwcTkSawOG2vcLWa6QFVqm8b\nBmcYE1qa+XNhpvdNb+GIxBy2IbtCl6Z95KhxbZbBaPdi5D4b3bGxC7ekd5x38FNmus7CsWj3JyJF\nYHnaXuFqJQjuuCDFH09CMq3gSBDnnHszOCYz/eYWjmM5JtyawureuDubpiizvTnmmmkM2eqBuuOd\nC7dYbcPki/iwvwtlmszFSmjgN9zbAYFj08r98XNgXt6eiASBL9D2CldbqvxhC149b5CYQCprlvtx\nl8zc2DGM4oShojfeGXwKN4bMg2t8W7ntIiunzzgtaErWT2QviNHdh6boVQ0tMXzfI3LlBr/HbEtK\nYnlb7n1H08r98XX2ejwR23b7p2MCX6DN75qW5SpOlVWU27xW4OSb5R5sjLDM0qcnHR7b+BK2lGI4\ne1p3z57lXTVEjLgTxI1ENovhbn/fhtATDRbSsn5AM6p72Alkk/1OI/ifCPuaB4s+FaYkgaw4beuu\naVmuUmbaLGb7ZUIXYvwUUGTWLCUhiVZSC7Ppo5zjhQfqK+JmSnGgA27rwO6tnVIPhFeDbKuPzLQx\n+5XmLNf7fmemZipdaJol696kGKO6+9CUVS42pfi+a8pbmzYPzJQKU5BAk3aibOVoC2F8Ya5SZspq\ntB9pTqEGsM8GKEIz7So5arAwodGYxlma7Kvq2lm21NsiA/tCcKBDN3ZtU7leHqipOqG7sEF5lKg2\nwrhqbMf9fo79ruy60xorvi0r0d2HppU7j6P4Xr9zb50y04P01YMfTMRDAl+h7SWuVjONUKXmWs9o\nWGCMQo7d3JvfAENaAZvCEAigTXawtGrIt/nIB8V6v6+IfWJlxYEOePiCga9Af+yB8EuYOLfrEUIc\nbPuZLcsFHAEX+BfK3MINhKuwDz527pz8tqxBdy7aZcCUZqjdskHxYS3IWFlQur2ZAuDSCvMtsehm\nIHZyJQhDJ+IhgW/Q9hJXp1Qp6VdfszQ3xlnOS4umax3IBETbH+CMYXPIsxu7EhzsYECMva9Au6A3\neYgSGwJDaGdRsON/W9ZV990gFD/fW+/NlIKen0KtwrxGIFXrAW1UL5N6h7kowtU5VU4Arjgs1fBd\n90wphkkFOD94l/612FZ684Qd0FSbFjOzFchUmg/CdkoeSUdwTeho64vJSoch+CzfWzdA2UE5PdM/\n6FeWQKJWOdoIaBGu8qmawFPNYEVAEuQvsNLLyYSsxiceEBqtxg7IdD/C6VQAgl+dwcEs1viy7o76\nmuSV7u4Wk5WMSPGfeOsV8vQXUmvTzV1RAqlaxWijoEW4yqcKU3gTTNbwi0/wfJjxYHy2Jj0bLg8f\nTAcBpjrWM6ZiDkjVuqs1sVLYHDrEiL3oyO5uKVnpOBS/jNs5/b0pHX69L0sgVasYbRS0CFc3qVpJ\ne/EOXe6Fcry6XgBJNCX4RdzO+S5BQpisV1kEErWyQLMaEdASXL1PVZZW8UaXHXUc5oXaIm4n/g8+\nJaX9CQJLcPU+VU9o/+3/KHyiGfT9hnLfGOMhDTndf12NH/g38xwab7WhYfgtgJxO/waBhqo/1KhR\ndFso74oAAAAASUVORK5CYII=\n", "text/latex": [ "$$- 3 u_{star} = u_{max} \\left(- 2 A \\rho_{star} - 3 B \\rho_{star}^{2} + 1\\right) - 3 u_{max} \\left(- A \\rho_{star} - B \\rho_{star}^{2} + 1\\right)$$" ], "text/plain": [ " ⎛ 2 ⎞ ⎛ \n", "-3⋅u_star = u_max⋅⎝-2⋅A⋅ρ_star - 3⋅B⋅ρ_star + 1⎠ - 3⋅u_max⋅⎝-A⋅ρ_star - B⋅ρ_s\n", "\n", " 2 ⎞\n", "tar + 1⎠" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we print the result of applying `expand()` to it." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAAUCAMAAABlEAJMAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAHr8+Dm6vLs9O\nn+9eft+PMnn1HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAxVJREFUWAnlV9mSpDAMy8WRkJD8/9eO\n5ADN1NIwHD1bW+sHoIHIsmIbt1K/aLr5RWcfcOUaYxr3AWClWvMR2B+AuvRAUK7tlGv7T2hjy9+S\nxqVOaTN0PxBx55WmIOt9CTuvXH2UhuHq0pvrggaAK+M9GDt4pT6yv40z5R63y6tLoTam3Ewb+o9M\nHaUfQCJaNR3Bbf4BcPuTmn2IQjvSmSn2flCG9eRjSDw9VAXAyuRWLSdrjpV/lsKISO4h2tSzXbpG\nZZ7DiQoNZWXf25WFKmmRxrChHWp+hcIk/cbJohRuI7rQdypqNTI6OWx4Onkr4/1UYl2VepztItQ7\nqGcp9IjmFKI2K2OzoqGZO1256zmc+kRd/KynkFIKhRUK+IESxW/SbOG+oSAY7w6b8fDlwH0+h2hX\ntvgzFERKKUpzXx7oE+W1LAKlRC/NJE3FDGvkN7hbFFawW5ercJbOhvfiNIlfQJy86IG5jmaO/ZVS\nApTkkrPoms4mI18WN3nFDzunGpZ5pMZiGAIWYzlBoGnmE3pqlAFwH3dNAb74WXOzw5kDz0Jq8bZx\nEZmnGlSvI3ZFpBmZ7LLHgJLvVGNjAjWTyCJmGxhY4/ts+5c276RhD6ZN0gi92mr8Pu6KAnz56Puu\nmRzOHLrGm+SFdXWydeykyfH/zw3EloFamYYHaNIBitJgMMBggjMfe8iSULo66sGpcDig6BaLaFV3\nLei9sD3AfVGAL05sLHRDoWcOZDdqXdNSnGwcOiO5TP1uILrcpsz1yA1jo8+hoRq59LxXW8IIfkEi\ni+xtB6ZNqdmSe1yg5DFQRt/K+iPcNQX6Ev+yPS8OE6k9FkOdKZjotxCdnStXyrqjMtpp5DKAg3JO\nsX9ip3C/CrTHauMZ29dUYYe4Lwrii/J0PbvfiwNJ1X644evPW08jjgglgxCkQMGiiuxY64z6nDVp\nNXXRGVwmC7ciRAxuCwc9KAztvHPeHkFkf5UKb/H3HtvE1scv1kH322S7npRO4IovA3lim/klmDi4\nNvk2XVLmecTNeE/c9Ot55sS6/+FVpNu/YF+xGSCTh81AmwAAAABJRU5ErkJggg==\n", "text/latex": [ "$$- 3 u_{star} = A \\rho_{star} u_{max} - 2 u_{max}$$" ], "text/plain": [ "-3⋅u_star = A⋅ρ_star⋅u_max - 2⋅u_max" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq4.expand()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's very similar to our previous result, except without the bracketed factor. Simplifying an expression can be accomplished by expanding or factoring terms, or a combination of the two. Whether to `simplify()` or `expand()` depends on the situation and you just have to experiment.\n", "\n", "We now have an idea of what to do with our three equations. We'll get expressions in terms of $A$ from `eq4` and `eq1` and substitute them back into `eq2`. For that, we can use the SymPy functions [`solve()`](http://docs.sympy.org/dev/modules/solvers/solvers.html#sympy.solvers.solvers.solve) and `subs()`, respectively. The arguments to `solve()` are the equality that needs to be solved, and the symbol to solve for.\n", "\n", "**Note**: `sympy.solve()` always returns results in a *list*, since you often end up with multiple solutions for a given variable. These linear equations will only return one solution, so we can skip a step and ask right away for the `[0]`-th element of the list." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAANoAAAAuCAMAAAB5wPfmAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMADp8eft8+v+/P\nXm6vjy5Oq/RRvQAAAAlwSFlzAAAOxAAADsQBlSsOGwAAA65JREFUaAXtWYu2qiAQ5SUoCvL/X3tn\ngxh10kRdZXfFWqeSx57Z82roMLY6uFhd/tpFqUQTvlb7VcW5luo/pUa8f9RWnX/RxZ/XLuqYVbV+\nXls1z0UXf167qGNW1fp5DeYZV2303kX/UpwxbeiMfbkPG9xruE04p2yy23TeJkslp3FlTM+3Hdm0\ni/e7AJ3ehL5lkzdxF1eeSdOe50A+eMaHrtpYvP7IEs0hOU1I2sBDs7Ster4PPWVxqL81Ojp3ypBt\nggkB3Ew4zW26JZvpkGKiRlV/lnn7IYkdGoSOCRTp8jR6zMJ1tYBnZcWgSos25MPRCsyNkzvL5erP\nBvFYC2hOisiuLLaajMx75hBF4nBcaNUBqBpQuGoTPj2ACJxHR0a2kjWwdXyZV/Z94KLz9YDKsHBo\nTMqW1AQoScp9JFwo3fmCmTTFQDnKg0ouXwBcrvFE7ZRRULNTjMdQtLFibhWhi3F3xsBCzwDlcryf\nRa2dAzK2OJJsHkOR1IH5NWeS/ugXsqQxn7bjneaXh2w7LJpAJakEJCgqwFwrE49nPHrIIpiaivYy\n+LYVkyukjwGIVgvaQB1668dutGPne91FbtZpgXDx/WjUGJVfEONDXG0QFQUgG3ttFZnLKFDLeCTH\nTSLIHFmlBeit024qR94oDOjTErQndUgDK/EdjqAycNdItBQSknY0Uq5WsiE6PXYjN0A6RmEAIzZY\nzngkp+VMTFHQlTcRuf+bIJu+TTUJLrFG29GJHsItEYlpEVVpkDjRu2mSNiwO7gblQOoe0IUOcynV\nbniQMw0eykAfDtSUqHOGTe8xs3yMQBCBWN8h9VBayFvYFTnfH3t84jpn4w1QcknxHY9zXuBNBsvI\nM9SeTm0+rKJh58eHDyAOscLSNy+jqNENArVl1LFgpnY0FNcuYVJaZ7wpPhNYDP2Mq9ojTZEp/Z8h\np3eO3DOkih0csl6gfOCyosZB7WDGRsrnmLUDwcx4LMpJMnWZwT03R/514W9h/kDsE498KCxNBSdT\n29ez+9WQfDM/VTBDIXapp6htsd+sdbU4TaGrIrXqFrta1psPIOtU7GXP7NnfTOKpOCWo5ojYyiy0\n2E+PfcGkVOi5+0jteYv9BSSeq5i+BOT068pdi/38wNfMooZgTNTKFpuanH3XkYT46VeZbzbpBlH0\n7AeuI58mBfnShOQt19EHtD23nv3IdeQK3P7qcGuxU2uOW3r1deQv7MVmjlxHLkblUZ2zryOP+J97\nPv068k4q/wC8ZB4pCAi51gAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\frac{1}{A u_{max}} \\left(2 u_{max} - 3 u_{star}\\right)$$" ], "text/plain": [ "2⋅u_max - 3⋅u_star\n", "──────────────────\n", " A⋅u_max " ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rho_sol = sympy.solve(eq4,rho_star)[0]\n", "rho_sol" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALEAAAAwCAMAAABpG2uKAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMADp8eft8+v+/P\nXo8ubk6v8XVtUgAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAytJREFUaAXtWO1yrCAMlQ9FEZX3f9qb\nIGpEzEC7rds7y48uJCGchBDPtGmyQ8is+F2FSsvWvyu4LC7RKf23EEMYH8TZu3yp8JPjl6Yz6+yT\n42xaXir85Pil6cw6+y9zPGdD/VGhu/VuTO8HY2/1qBjvtzeNmti9X1VaHhLvVrMpXgy/u1a7Ucmx\nq9252zsWUudZ9e6lbEKopBhE2Z6r1cKmWPf9dUu5ZDrXFKWS41lV7lOxiCZhvkVWtU6R7K3Ltamq\ncD0tjKGyDUWsutqbZBA3PffgGVDLJQvEGHSj35/IqDtTeQqH2HyxLAamzXQAVu+I8YSZrSESa5xy\niOV4tS+R7IAyxuhS+xiSHmDVceYZDxxibRpfMzb3EYIyZKhVqaXWWvq1bESPyG0xYh38tW34IZW3\nv7wGEH9pbBA6MlZHSqNoioitxzhk+FtykAr+pAw/MQW47/uI+/1hpTjWMlPxEyJDM2oxL9AxsGtA\nc123iOgCf5NuwlYF16VSMGRtyI0RMZwdcUTELf6jJpTxNA+znQc3dUOAbMdOYiBumo2esdqPwSFO\nT+70RO7i8JHOxpsXq7YM+IBBeUSMjUVZ5aHH4YM0GNQMaDUqIfRWqbM7DvFw/tgC3K4taZ5JUmJA\nyvg1t+MAE2hr1k92Xta2YQGfwiJpMSctCOWqCMLoYf1JERMqKfzpc+Q8BKCLvoPh2NMxuQWUsYpl\nEvAhajc0ABnfIuQWN4VQTrtTxESJxmQoA4jnIkrAOCUOQxnHNcaIWZVWQPZ7AbcJJQEEReoGJWS4\n+2sOJUVMcSq3UkzkydKcbidRxqXaviOwFljYBlDbZYStVuKLc41Y9LzoM+C8syDtzgUfZK6Qgbrz\n9WRPmYv7cHb7VSiWI08bwXJEdt1BJe6uwRGj+9slRhVTfQDeCJaDrK8vuMLPE6YbwXL4ecyUyhOY\n2DM3gqV6ZEV/AHEtwWKj/xVlLcH6FVDsIYRgsXbvozwI1vtgYpEcBIs1eyMlJVgVtPvBCAjBqqHd\nDyI+CFYV7X4OMSVYK1cto93PIT4RLCSwhbT7OcQNJVgVtPtBxOTo19Bu4vAnpv8Ak6QYiNYLWxoA\nAAAASUVORK5CYII=\n", "text/latex": [ "$$\\frac{1}{\\rho_{max}^{2}} \\left(- A \\rho_{max} + 1\\right)$$" ], "text/plain": [ "-A⋅ρ_max + 1\n", "────────────\n", " 2 \n", " ρ_max " ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B_sol = sympy.solve(eq1,B)[0]\n", "B_sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The arguments to the SymPy function [`subs()`](http://docs.sympy.org/dev/modules/core.html#sympy.core.basic.Basic.subs) are given as (old, new) pairs, where the \"new\" expression substitutes the \"old\" one. So here, in `eq2`, we substitute `rho_sol` in place of `rho_s` and we substitute `B_sol` in place of `B`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAAA/CAMAAADT9iBqAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAADBxJREFUeAHtXYm2tSoIbrRpN7z/015w\npNLSsv5ON9c6pzIEJFLUL3eS3JIqdouYT8i7LVDld9evb4u7RX7y3miBMr23VsWvv1fgJy2WBSqW\nd7F4ReBT1FUELv4shi9W8DfWoyjBbav6Sc1O9buz987rRz2NTxlvC/TTmCTsUY+vKb21P03Y/e5t\n4U8r/DFQFugG8NxxUpdPOBY/UOmm1Aw3CfrEXGKBO1s5jwrc1wVU09fkejwQILmvMdnVh4S2T5sX\nKn53jZqGdtdOHwFaICPu8q8tkurZp768c0jkU+/spkFaN2kj+KhlpSkaa/bfzizyYciJVzDe5BZs\nnnu6jksxngwz2VH2WZKcf4CeQv3IYniUj6RmIk/Hp8CSpmNN/agxwlLBY9dF2SdFaVZoej4aKFif\ndEPE6e+lGF9lC6FZ31RVBc77qNTeM9nxOzuLUVQde6Hn5hOsZI6T7k1K3uQ2OO1fTPEezVKMtwtm\nuNDa/SZIT/PcfLojrhojBAvJGz23wtmdalLzLt2P+9Q0oesO8R7NQoy34yZ9vLfHX6gfZXfLu3Q6\nWMDKvNFz+UNKseHlKRddU1ljbDXgdEwXr2ERYgIZRgxZZB2jHW4JF34xZhZe67mDCRboVE8NDfCY\nNpg1irb43EPnYkIZDuqlOif7itJZvD7JqV4fpWF/p+dWrFWxQpK0ZABfQUtc5EmGd5vTnbYUE8yw\neVp4a5xs1F2VyYt9xnDh+3R6p+fCUKxpVUhAl2taaInTLqmxQeb/ThqQiwlmyMxrdVJ+9OKFHh5E\nZ60ZlnzEoS8PnrzVc3EWQc4ZEs9t0GMBo8WD3ZABbjeQRLGJKMbB0EzLLZ/Ngz03wXDq4hRnJva1\nngtjMRklGM9NZXzJ44Q06M2HuVedZg+Wi7Ex7NzByJM9dwgyy8wUnheRmvUXem4nhq7DJEdmGk/H\n1107aDF5nADeho1nVSQd/MHUtjB8IVe48Aj57kTFUIbACuKUomIDL674wYUSkbCzE/FurU7fiROE\nbqlRTVGq/0LP7Sc+6VKrtlaB73veBuOqMPdp8DZw7Xxsx3Rs+7xqueumWcUBeH0+Dmzcmr2hYgjD\nZMyrlMHbMDD0XMUP5GRSBEzN0dmOrYf8D+7l6oW/THYaR4KP58YYCcaygxp2bfEreVuq19AyMZLv\nB4YJnREhUT3YDxws7XAaCAOLARvbEcZODINhoKi7bnMOgIgxDKEYNOL4jtSoheIHcgDM0sg2vHVZ\n1JUPnK5KS4NW2m5XSWQxpi+G4Te1g4wHXZo+CWiVGKSVS13oqbOSZQawJ1tOvto68aFzOlTpmDX8\nM7AU/JSHpNzTarBFI+yxEacK0UQMZZhNLTaporjhh3JkKlx4k01Ld9dMAi8N2svJhTTLKGhJ6R7h\n2KjOUPC6DvMlgFYoJY/cyaXlwLYiSauVFNLKelNlFhUNUblPqlv8yKPanocH6KfoVX2LYS8OT6Ct\nRSLu0rNSqwsjxjDsig6CD168KAg/+T4ozitekGEsbbtbRp5KUx6zNKhAdpTwrCs1O2NT53geX8aU\nxa/EfAmgFQrqTDN2XG1TshiyLgkfZUukleGzf8a23jh0a/SqJoUFhQS69KrGKOKXwDob5oQm/IIW\nKoY8oc1S/GTwIJjxuGTN11h6fY/iMGx3Q/OIxywNOiEGK+UTI3wSMZT1Lj313CsxXwJohepkcSJr\nVbMBR5jcTCrH78iRVn6kimpwt+wFxr0DeFpaZjiganBkhiBJNpbsgOMmI8TSPGIugY3ml3A5Qh8X\ntNFYWulNj+x3bqI1n1WGeszCoPyRtHz4n7vCGqpY8Hk75+oz0gqWAQUk0ArOGI5p4qWUL5CP4WOS\nA0ir3sSY8SpwlFPh+BLCWNrGOS+Gc3jUVcejPWZhUO5Y4lHvfi0WiDQSFVu0VloPW7VP5EmgFcSB\nTEwmxQJataYJCaz/AaQVKP+Y5ArttaVtmsJsBfVcnIEOTG7PTeYGxc68Up67bbdQpJFQ+SbPxUid\nJ2i2eHWOaSuZmAMspBSMYUAYDN16MNLK1C/8TFvaVhSeQmaG5BmrhuVslq0Uzdvw3LlBrZ7bIAxe\nJdWDBSONhD43ea4CWuVgKfTcg9pSG/Lzamoh8OrhdQ/m+GCk1aqaARnK0rYiFUTOTHsuelowRHPD\nc+cGtXquTamj0KWbPFeaq8f4Hj03GBdlrTMOlfn6OMC9gzmGz0c4dHhWtnZMi1q4IsIUDoMvpuyG\noEsuG547N6i/5zqQRkvJy+sNz7UDm1Rb73mU8qQ9+WISeu4hbS0KQZuLEuBxODh6IK0sbJGnZwUf\nQ4Y6Q5KWtlWKNTBh0fBgDXq9H66UpFuOzvnpf/DJM6S65gd4hDKZkdHaczshq1IviyqzPNqQRkua\n1fWG5yYa1QQnq4JhGcJAOZqOTQPDiU+c6gsDWlkUkoAhhiaycdxYwTKGjljPMKtcQa1ccV2pjmGe\nAhQI0zf+kK6Oc2zwU+NKg36w0VC1MAbFHGxzkx+PZcX8D2biXJ9OZj6IIo1g0OgHXdryXJQVKQmg\nlaj81CAK6pC2a21q2eZCGEI5/nGk1bqe3jka0rYqIdATnVyW5e95UuOKmrezAO1WtDDDbnHPFdt+\nmd2aHJ7Lm+ZQ6NJNnkthTaILOaTt6mnIVz7DpoNw/OtIq3U9vXOopWeFVL8pPZe/5zzMDcK5bXju\nXDL33A5nyoratK4zjfSFQRoFQJd+urHnfEzbr9lGOZFAK85LeO4hbde6FD9obYsa4y7DEdrfE0ir\ntRCvHNsOzLY8J7MgYicXamlK1KkWUYwMOo7lwomIAGcBfhueO4euie/6qyEdSxyWbyaKNPKGLs02\nDvDCfG3q4LppIKpDO/0QVnZIWwv7vm0Yd9wZxxNIK4uMzSwJvrLtwGzLc/IKInZySYylKU03CHhb\nkrVwAq4En8mDT+EgDc4hGvXFuS08l3jMArqmmsQZaImqRM8N0oiDP7ygS7xRp0wuOr8MaAWzuWol\nyNT/BNIquP4CfGXbgdmW52QfROzkAjfWlrZRQ5jbqWF3CM5t4bmENbobSYswlNzZPvWGLt3lue4a\nQ0W8td2utLl7HGlleHieyU1wbDsw2/KcXIOInVzgxqaldUEe5sqrEPP3zkW3BXTtoOf6Q5eu/15I\nmudOoNUJpJV+tL4nFHxl24HZlufkHUTs4rJhaV2E7rbo7yy6uOVkAV2L9H2jRZDK4jOh6uLK46OA\nVgA/DMaa2I1DwVe2HZhteXZOkBtE7OYy77atdKP/PK61/CpzaVCzH9uKNFJGFb7FTSAmSyn6JKBV\n4kJaKWV9jxR8RXdgVigsmmfjObPlHrGNgS3Px9LOTt/G0CNvadBRzCB5lDxKor4X8i8fCeXlL/DJ\nlDAbp8BXdAdmhcKiebZqzGy5R2xj8Ni8KN83btcudCfYYEzWtvi/fZeAr+gOzAqFRfNsFZ3Zco/Y\nxuC5eYvvG69QlCKNffgHY7J8mP5VGgO+ojswKxQWzbPWkNpyl9jK4bGZoW51oCKqs/Mt6sBk+RZ/\nFd0MfKVr5o/CerEtxUqdtskVJ2ZrY1/uNkyWb9lX0c3BV7pqISist9qyv3yjEFjzD96eiWKyQmBG\n+tG+5WQOvtK1Iigsnec6obZ00fzF/DzK3rY7NSffIO5QytsEkxUEM/Lj/neoFuArrbhBYeks5wmx\npZPmL94o55+UX1OF4N9RMZisMJjRNer/M64L8JXWw6CwPJCvxpYexFrE809O/96TTxWDf0clFsrL\nR7nH0qzAV1pTgsLa75KILfeJtYjHn0T5vaf9Wgb/jorBZPlj0vbVeAuFQWH5dEnalj7Ef8ZEUX7v\nab+25jOhfdolRQjMaFn2pdcUhRWCfA2CyT7cduL7s+uVPC4nDszo+hreKIGisMK6pBCY7I0VOiDq\nTFsYJI7d9NPYQUr9VeIZCiuoSwoifrR5hsg7nTorKz7jct7+bgRZgKCwgrqkIOIgje4mrm74GT9Z\np/RrdO9+um+WZ3638/paDmI56HpBn4T3W2C8sx3k38W/36ZfDW+wwP72ClGVSPlWMVFZfsz+nxbI\nPD4mimmZ+RaSMTl/vP5XFhjdexBeZAecTfzSZ4GTFuh9vjw+KWNRvCCzOYtb3+VnAV8LwO9u6/Qf\nXQ6BiZm7U5cAAAAASUVORK5CYII=\n", "text/latex": [ "$$0 = u_{max} \\left(1 - \\frac{1}{u_{max}} \\left(4 u_{max} - 6 u_{star}\\right) - \\frac{3 \\left(2 u_{max} - 3 u_{star}\\right)^{2}}{A^{2} \\rho_{max}^{2} u_{max}^{2}} \\left(- A \\rho_{max} + 1\\right)\\right)$$" ], "text/plain": [ " ⎛ 2 \n", " ⎜ 2⋅(2⋅u_max - 3⋅u_star) 3⋅(2⋅u_max - 3⋅u_star) ⋅(-A⋅ρ_max + 1)\n", "0 = u_max⋅⎜1 - ────────────────────── - ──────────────────────────────────────\n", " ⎜ u_max 2 2 2 \n", " ⎝ A ⋅ρ_max ⋅u_max \n", "\n", "⎞\n", "⎟\n", "⎟\n", "⎟\n", "⎠" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "quadA = eq2.subs([(rho_star, rho_sol), (B,B_sol)])\n", "quadA" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAA2CAMAAADDG3zTAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAACJRJREFUeAHtXIuWoyAMRUV81cf/f+0m\nUTBYrOC0FdxyzkxphCQ3pgEhIsSvXG2BShXd1Tq8T/690LzPLhdwAreq6v4CwR8ReS80HzHRt5j2\n0yiEqr8l7sNy7oXmw8b6MPtuAM8apw9L+Rb7e6EJsRrcxVgKHwCb8pRW90JzygSxdMr57bxYqSwz\nCvStNPWAyr3QBAC/vqkshqFgN03Rj1wqm/p3Pbdy/Djm1dKuL5mO+323UtJGs48zgSuy7IUs13jQ\nD6i0VL3ohscbo9dWjqdp5KJZnwuxxq/dzlspaaPZhZnEhWIqcHLcaGVLClkNLh7J6Y1PY1s5Wt7R\nZw7qCdE3VVWBcx2VrZS00Ryhjft69QBXqiaKVKBp9yB1pwlda5jeF7Q2cryN0qN3d48JiodnbaQk\njsbbSPE2zDBwUSnmB7CyxlnNMMEsp3ufd4lZThjH8CH5XmiWG5Pox7AOhopBqCGAjVmDpHGOZezi\nmSrJCeQ4aKf3FngvNN6w42tYqVaPhUK0bJpcQSSThcjxavP3KdciJ5Rj4zEIMqPeCw0DlmRVNq0e\n8nD806WFSJZ1osaARv/0hbOfJCeUo1rd3k/uvdD4YY62FTwFLqtFzLMa9CjYCabJ1sRi2RGMbmDF\nylBAOTsc13WPDftgz8Jn2huh2Zgjta+D9pzVs7JlfkPjYEbPir6oYIXAFLsPyXFx7HZH23DPgieP\n5XdwCzS2BZP51j1a1HWYlpn7Q4+GtK/SQcShcRCcAYNPJUUHf7JaIpFcWuMn0F8ULodzBF4wEMtK\nDdRfM4QvWoYK2C7kUlCZtNG8MGcCl/qJPKvWv+5h8bCefvW460M+B84AF4qxHbOx7YuqJdfK8qrB\nWVBfjIMaidEeZC6HcRRjUWUK3HVQ6FmaIQjKFxlCa7THmdO5FKTrvmmi4cgSrJcUi8wafD4/ifWD\nwoLO8gCX6sEZ4P5nHS6d4lAzYLAawasUTsagRd11r5/hmJyVI/SDKIhOXKMamiEIekjRzEGwDUlW\nYFKAn0gcDUJIt8i8VDne67kskYfWvCdamc+GKhvzhpKFM/AjmhGRI9TgEg3FtpmoeTg/mRzOMZ9a\nlD1Ps1aGKGgucno5yupmunV6aLI85wkBG0A7X7cb7zvNgshZOagQWx8zlxWfIpHPWJ1oVtXT8Id+\nhDe9b3HahbN6iFXYmFzO6vX8ZZWzcuxkB6Mr9ZeSMVwcduH8zGufskrBNvGjKeGHVeln2X1Y9pXt\nxrt99dQ3OeSdOPG05C9MmfDl6IM3Cm96k8F6p4ARq6pxlHwIWKdHSnCpYVQFRMgUfriaIfcIGneD\n+eoO0aPJ6JmYFni0zh6f2413jy5PTQr7fg34oDT9IW13w+9JHkyh9iOixHnXAI6QlTm0yhqcuWMS\njhpLZSvqYOwijTCZoylbCXwMQ0GCqL1PgoOLsabFjqalJ98iaMiHIGenEWiwQZ/2jy6j5IMxZE67\nkWbz21ykr72Z4riufpkm/fL99rWKHc282kOr0fsg3FfObO8zTrYntOu+cFjCgOFo8zNkXulfjYe8\n4Rfqf59Sxo2m0p51wuZntvfZLbM8QU6DVArnJadTECx+TM6veokFnj2rwUw0XfaHj5Pb+wyk5QnV\n1MJspoeMpdCEAcPR4meov8pFFnj2LG9FTm3vM+6WJ1Rzqiek5YUmDBiOFj9D/VUusoCnZzn3+F9s\n779GAy/PQKlr+piHYYhZ2EdN4FjBKQgOfq8ViPqq09ZRa+xUrtPzrHnB2dkGiWaDHyqm0e72vmnh\nrnTErcFXC6pldxbmWdhWoTquhAGxm4oCvRz8gKpH9FQ+ET8VMsvyT9MSRPOgqdT8zE8wcB3GlJ1F\nAL7xbm3vG0N4VOzRq15iFsy2XBz3U1GMJJufIf8qF1mAdvT5kRU+nsU33inozQkDYakoticoWiPN\ncVuFczxORTF2s/kZ8q9ykQU6fIFE1jvBaU8ptvHOtvfDUlFsT5APiFayxlkX4+iRimJUtPkZ8q9y\nlQWqIRvL0P0L5/Z+YCrKxhP6tlHkWIInDBynohi7bfgZ+v9ZcZ2Q5qJ91DpWQoCvpHXjfd3en/MC\ncCPSJxXlyRN6nZjAOB6moqz6PvFbLx3WutCf1iHHqxosSFwnpLloV6kZKjcoFaWHYfioHKeirBx8\n+K2t7VoZ+haN3T2ibzMS1wlpLlpEir9W5YJUlNcKeV5dX9H37BBtswWJ64Q0Fy1aHBvFLkhF2Whw\n8qt6rPvhJ1lE0o0jcZ335qJFovod1Sjk8Ie8sJgswpG4zntz0WLS/2a6wHsQ3LPw5TCvcjLbx4v3\nuUYcCT/vTUPiNKeE+CA51UyFCGtouX51DKqqGjyeLc5n+3zQLAwJP+9NQ+I0pxaBJ544efyIxgK4\nEaqMZ+H5MF6n05zO9jGC315hSPh5bxoSpzllRwjJqWcqRHyZEHIsZnXpPUSvDNvT2T6fs8uKhJ/3\npiFxmluHCCG5FU2DqhrYf2/mlA8hH+hhmYlgLyCcyPZ5we0dlywkhmHSkAyKBCudwnyVYvGs+fiQ\nxu8QEWe2z3UmsJEYPVKGZECkWJlfrO+W00zJWUSNK/LwhHhwhAjP9jlu/XHj2EiMOAbJ0HYrkUHa\n1TOBCzqPcfEssixNszzyNli2j0frTxtjg8SIWyEZ0n4lLkj7esZ/pdPnDM3J0h0dL4LHTvrkbazZ\nPj6tP2yMDRIjbYXkE4ajgmQwpFfphvnkEJG3UIHlBniDElKK5sfE4yNEeLbPceuPmucJiZHGIHkE\n1oggGQS3qMCcpNNp/h55Gyzbx6P1NRZaIXkF1hQgXWPIv0mlOcnCIixvI6z137QM6s0hhQXWaCEF\n4Y+jcafXS0GdsLyNsNZfRMsh0Vk76F1eJzlFC+mL1nubqNFvHett8r7AyIL0i0JfsLhbhNdOtLtr\nrFQG6atR6B9wnYKbM4kVAAAAAABJRU5ErkJggg==\n", "text/latex": [ "$$0 = - 3 u_{max} + 6 u_{star} + \\frac{3 \\left(2 u_{max} - 3 u_{star}\\right)^{2}}{A \\rho_{max} u_{max}} - \\frac{3 \\left(2 u_{max} - 3 u_{star}\\right)^{2}}{A^{2} \\rho_{max}^{2} u_{max}}$$" ], "text/plain": [ " 2 2\n", " 3⋅(2⋅u_max - 3⋅u_star) 3⋅(2⋅u_max - 3⋅u_star) \n", "0 = -3⋅u_max + 6⋅u_star + ─────────────────────── - ───────────────────────\n", " A⋅ρ_max⋅u_max 2 2 \n", " A ⋅ρ_max ⋅u_max " ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "quadA.simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's a little bit ugly, but that is quadratic in $A$, and that means we can solve for the roots of the equation. SymPy's `solve()` function in this case will return a list with the two roots. They are long expressions, so let's print each root separately." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A_sol = sympy.solve(quadA, A)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8IAAAAyCAMAAACgT5IbAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMADp8eft8+v+/P\nXm4ur49OcUBLqwAAAAlwSFlzAAAOxAAADsQBlSsOGwAAC39JREFUeAHtXeHCqiAMBSUtM/X9n/Zu\nQwQSlaklfRd/lJqMcWAwBseEyMeIQNdnKDICP4JA2/2Iot9Usyi+mVvOKyNwCIHcXGfwNa/ZrXwj\nI5AuAq/met1kSsNedZPXI3K1BvWQj19BQMhbdW17KVVxG65Vwc1d3i/Gw1XmqvPyflXOOd8dCFT3\na0cdWZUqIRN+peQR7KjOU5I8Mwin4PgtIcXlc7+ETLgaym/hnnA+jxyST7h25qqVw9WuY0ImXOfx\nRwiZkFc0b6/5zhyBop7f++qddEw4D8JY8Tkm/9Xmf0Jmlw/D6ZhwfXVvdkJ1HheRwirF8VL8VxKu\nbrjJmHA5JLDEdnLTq9STOb2Xw7UBzpMBOCKOD96R3A6kbS6O4SRjws+/13jBfqtby2oc2Y82cO0A\nzyT98rccnl/O0c8uGRN+/Dk/uh0guKxuPt4bV9mPHgHaA94Gth/7uX58THSM4FRMWA4qRt1feqas\nwYR7VoQ5+9GmgneAZ5J+/Vtd60CmYsL95ctrn6l53sp//+dckUOo8sA7lNWRxBW6W9cdqZhwd7Qn\nuxRFv/6c6W/7YIWnis1Z1UeCPAlhJ5bBS0hLR0moejlcSDqs6/vwqFMIBdcHtwZ3Pqq+TX35qpkA\nbV8sCxb3rQD2R4I8KWEnFsFLSUurJDWte3aeAIfhGAqKumip6vrJM5oN65bPPRK7cc8dUsInc97I\nCX/ut+IiHwnypIWdWABPa4koPU+NmjR7KnhSEvURomYFPHSaP/cph0O7K1vqAKRqRVnfTxyP5asV\n8sVzhqFupE7RFlVVcXysTT/6E0GexLBbAE9rie2+vJ9pwkUhhWJXsFFytMPi6CxwlPPTX9WxgPSL\nBuEC3VA58JZxVmF74pJfz+9eOpzTlnekvXJMeNOPJmVPDvKkhp0Igqe1xPJ3x5oKQTh96BDMjVNJ\nOi0pacSoPxqLNeWL+j5mwiPHdqBtMvVw3jBc3aFvqPhOfrurG9n0owlKZoRsC/7ksBMh8CyLWrUn\nmrDUVfviB2I8JbMJQyt7vvdjZcWY0z41Z/NFL/2oUVR5nhmLhnbf8CQyvHlbzlU/2gCyGSHjKSrS\nw04EwBu1FKJVQpsws5jhnmzsnckEmQJdJauL92eFC/flu+/9WKeqOt4KX+7s6AZdat8UeKvn966B\nchMJkimx3lwdmjJ6TFHotc2YBpDNCBlTUZEediIA3qQlhEzIhLnFnOD2TiYTVuwm4yp5zIf0NPrd\nizcTRnwY9vdwwr7YI8qn6DDAVezyZz0UK0U7P7kSi+jZVTHROzzfzFMCYp4jIJsRMq6iIj3sRAA8\no+UT+nU0YXYx3/A0l3daA3gNii3QVTKbMODpm7BCYKt339rAPv92n3xAP92U4oYRbvqYP867I4tH\ny5aoYtfIVDWYN7d0y2ZvANmOkLGLnh52IgDeqGWLzg2aMLuYC3XeY/hE3gfFFugqiSb8K2/sO1lP\ni6tnwvKOo2rjNi77ZOjMeZLe7AsMIbzF4i+WtXNMvi3mBjFuuSBxcTXCreBR43AGrXgYb3/yo2dP\nMgDhKirSw27FhKmTQxNeKGaocUz3ZrDiL03d90814JaZUJNZrGBPyTwKA5KeCWsCZhFPw7TNsBnn\noORC84icsIY7HVO100mNvm5IYrnoqAdMWEzS4cTJwBTd8aPfn2QBwlPUMeFUsPOsYwRK1/CzUHAM\ntVqoDwfV0Ok7rPqZqoV1KvDPmbi5NZxNGKA07ZhQJTDFDX1RiEtjKBZetanxlmPjx28bs55eXksb\n30p4llxokIPJ4mRo+f5nqadKNfb7rkTQB+pcVqomHYxScGH0VMY79gWGrtpx6XjFj3YACUnw7/EU\nFelhJwLgaS1LMsKhQPTdYu6vYQ1djV2xK3C7gj0lswkDft6iEoFJfs2zf/RN/2iflY7bNl1Fbxtr\nn32t7DJqDSaGR0thLdxiiUaHlQJfkTJIwNtHO1C044ZjgCNR9M+qUdC31ApN2CgFGXWjnsJo9CYw\neKn7CTH50fOHLCDz32Z3eIpOmqaD3aSSUzQXTyqgWx8HarjD3TSSIoo7cUMlVxaVPkJLgSwv53xA\nL+ofbj9W0m4ojEGWTYkuDiJc4/CLXDyFcSocFctyCv+MA1hbo59FoR/chId7AMDQImX4+oxXLxrF\nSR8rEfIG3wB7ixv+bJSCjOC14LBfDw/Oq2T1jMFuXiAB7ocFJMal4Ckq0sMuBJ7romgTtsU8UsMv\nNGH9/nIrMKaCPSU9H9KtOfERWgrkcD3n443pgf0YDpv6gK0UTT++vqIBg6UZJ1nLDeymQNsZb+rn\nwYo0PYA2NA604aapq6bvCkAQRskoGUaW9y27l+qwbjHuYSV2wwPv6amwVQoz0geLut/QILDyCngH\nkIgBh6VogtjBoGh3uxhAHY/rMdyRW+cWc38Nty+lbrrpuQI3K9hXUptwgGTzEVoKYGI5HwagA98B\nvWOkGTqKedalOeAUdpzyksFiBbUPnNZigAtGX0xEtWZS6+FwusITmj239CgafYQML/l0ISua7noS\nS1mCc09CpXSUGnsXeBaziz4kLSutjNsWkKgBxxZ9W1FQkjpHT1krIBZ/L/l0sQ+7MHgf07I19es0\nmQjcvBommkOIZPMJWgrAazkfE9b7T0J6x0gbuTz2UWcfMs38xl+w5tAyigbW3gU4qtUNneO7gN1X\neGc8lB3DzS37HSnDJtg6w7/C6rRiMO82SrnGQG7/lpjpd1xWWvGjdaBlfJq6Loq+YPdkPYDRJZhk\n0sm2ovBYatjpOZNfkC9rGYGbV8NENlwk2ZxMSwFkLOfjHaYd14t6b8nymB7wsKX8u6u5Ej3kGhpr\n8+owcFRgGAv5f6p/KWvB8Mzc9zIaRMswCTa/e5hw09T8BbpMSgEBzaRkMQyheUK4bMWPdgHhDYub\nipLCaWEnFsD7ppbbuPlKEuV/iWRzMi0FqmytuzdNMP57Se9NCc4qKD1rX7xDG2Y207890HI817e0\np19K5rs6cFlpxY/2ADndpQC3LCXsoHsO98ZJaekrqV+845BssEV9ipYCog3ng0nMWGjnjt5MgS7T\nA4Q7r7+bhasX8vZuA4MlmUOFG+GyfvfHasfqAHK+SwFapYSdWAQvJS19Jd3X3yHJBo+P0VJA9sj5\nOIfpQdrClAz05gp0mR4g5g++hHYEJ+KrGAK7GSLS5UeSQMB5Ca1ZIf4cLQVKrDkfbGLGGla7yEEu\n0wOF/71Xwa9B5v/WDLQr3L+Zr34FAedV8EiygeODtBSQrnebsokZpNnCxy5ykLvHFOX+wT9kWYBr\nflsOgaXQ+WP5TpIIOH/IQiQbcCnjeToL9IrVgmoT3pMyyPSAvNbIQXFMD1TYi7uuluAP/viI31L9\nB0v/40WybBpDFDlMS1lHRJswGB4uLtrc1xPpX8NMD6N3SOAylWdOR7n6Px5jEPjUM7hvPx8/isDU\ncCeSDZmC5ulsl4m2QeDuHXg0jq9hmCn8lGFtJr1dgUymB0nOfxEeBjjfTRuB6S/CLVGETCH2lRUO\nvSJi+yxiYTgf/JRBJK3ejsAIKs+kh5VK76iyl/ksI/ALCBB1DhS1JJu9tJSo7bOQk+F8WGJGbMoQ\nnlZvYQXCKtMmlSewlyEPwyGE8720EZgGYUuy2U9Lids+azgfLjEjLmUISqu3Rx1hMj1GyZr0Fcom\n38sIJIpAYAc0Tmw1TydqdOSzSgzng59yA0QrkMv0GAVLM0/fyCj/nBFIBYEKWC7vB02Fx5u80TFy\n+2yAmRKZ8l3VxWsu08MIquhd7uYqf2cEUkdAYlN/O7zlURbTNXr77IzzEZ3yTdfFSy7TYxLU5PXR\nCYt88gMIjG+m8DT9NC0FMkuA8+EzPRwAxu0tzp18mhFIFoFwc3UW+U8fHUcorud8+EwPt4a6y9/r\n5WqTzzMCKwjgv+OMxz9iIZ+Uzh/PegAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\frac{1}{2 \\rho_{max}^{2} u_{max} \\left(u_{max} - 2 u_{star}\\right)} \\left(\\rho_{max} \\left(2 u_{max} - 3 u_{star}\\right)^{2} - \\sqrt{- \\rho_{max}^{2} u_{star} \\left(2 u_{max} - 3 u_{star}\\right)^{2} \\left(4 u_{max} - 9 u_{star}\\right)}\\right)$$" ], "text/plain": [ " _____________________________________________\n", " 2 ╱ 2 2 \n", "ρ_max⋅(2⋅u_max - 3⋅u_star) - ╲╱ -ρ_max ⋅u_star⋅(2⋅u_max - 3⋅u_star) ⋅(4⋅u_ma\n", "──────────────────────────────────────────────────────────────────────────────\n", " 2 \n", " 2⋅ρ_max ⋅u_max⋅(u_max - 2⋅u_star) \n", "\n", "______________\n", " \n", "x - 9⋅u_star) \n", "──────────────\n", " \n", " " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_sol[0]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8IAAAAyCAMAAACgT5IbAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMADp8eft8+v+/P\nXm4ur49OcUBLqwAAAAlwSFlzAAAOxAAADsQBlSsOGwAAC6dJREFUeAHtXeHCqiAMBTUtNfX9n/Zu\nAwQUFNSSvos/Mk3GODAYg2OM5UMiMIwZiozAjyDQDz+i6DfVLIpv5pbzygicQiA31xV87Xt1K9/I\nCKSLwLu9Xzee0rBXP/j9iNytQTPl41cQYPxR39teyqp4TPeqYObOnzfjYSpz1/fyeVfOOd8DCNTP\ne0cdXpdVQib8TskjOFCdlyTpMgiX4PgtIcXtc7+ETLieym/hnnA+rxyST7h21qqV092uY0Im3OTx\nhzGekFe0bq/5zhqBolnf++qddEw4D8JY8Tkm/9Xmf0Fmtw/D6Zhwc3dvdkF1nheRwirF+VL8VxLu\nbrjJmHA5JbDEdnHTq6sucnrPp3sDnBcDcEZcPHhncjuRtr05hpOMCXd/r/GC/daPPqpxZD9awXUA\nPJX0y2c+dV/O0c4uGRN+/Tk/up8guFw9bLx3rrIfLQE6At4Oth/7uXl9THSI4FRMmE9ViLq/9EzZ\ngAmPURHm7EerCj4Ankr69XN1rwOZigmPty+vfabm41b+xz/nipxCNQ68U1mdSVyju3XfkYoJD2d7\nsltRtOvPmP72r6jwVLE7q/pIkCch7JgfvIS0NJSEqufTjaTDpnlOryaFUHBzcmvwYKNq29SXr9oZ\n0P4dZcHsuRfA/kiQJyXsmBe8lLTUSlLTembnCXCYzqFQURfNq6bp4oxmx7p5d0TiIPfcISV8Nued\nnPDncS8u8pEgT1rYMQ94QktEqbs0atIeqeBZSdSHsSYq4CHS/LlPPp3aXdlTB8CrnpXN88LxmL97\nxt9xzjDUDRcp+qKu6xgfa9eP/kSQJzHsPOAJLbHdl88rTbgoOKuiK1gpKe2wODsLlHJ++lSfC0i/\naRAu0A3lU9wyziZsHS75jfHdy4Bz2vKJtNeFCXdbs91dP5qUvTjIkxp2zAme0BLLP5xrKgTh/CFC\nMI9FJc0/+7+Qkurn6o/GYlX5gs7nTFhybCfaJtNM1w3D9RP6hjreye+93Ui1MYTs+tEEZWSEbA/+\n5LBjLvA0i7rqLzRhLqr2HR+IsZTMJgytrFv2Y2UdMaftBGfzTS/9aFBUeZ0Zs5Z238RJ9HrzaxPW\n5dz0oxUguxGyOEVZetgxB3hSS8b6igkTjiymuyeTvTOZYKRAU8n65v1Z7sJ9+e6yHxuqugm3wrc5\ntD2gSx3bAm+N8b2ro9xEgoyU2Pj85bUJv+Yo9NZmTAXIboQsUlGWHnbMAd6sJYRMyIRji+moV7g1\nm3AV3WRMJc/5kG7Vfu7uwoQRnwj7exlhX+wReccGDHAVXn82GKC6op2fsRIL3+xqZcLFTO+wfLOF\nggqQ3QhZrKIsPeyYAzylZQf9OppwdDEXeKrLJ60BvKcqWqCpZDZhwNM24QqBrZe+tYJ9fTaffEE/\n3ZbsgRFu+lg/HneHF68+WmLlWyNbmnBVT+rNLYPP7AEeCYgzQmYVJ7ro6WHHHOBJLXt0btCEo4tp\ngaQvRgyf8OdURQs0lUQT/pU39l2sp4bSMmH+xFG1NRuXftL1zXiS3uwLDCG8FcVfLBvjmH1bzA1i\n3Nwj0bsaYVaw1BjWrOF4POg0e/49eylvf/ajV6pEABKrKEsPuw0Tpk4OTdhTTFfjmO+tYMVf2mYc\nu2rCLTOuJuOtYEvJPAoDkpYJCwJmEU7D1M2wlXNQcqHjiJywhjsfc7XTlwZ9XZfE0uuoO0y4JOkF\nrhXXtdFHqKIbfjQ9Ij9QgyhA4hQ1TDgV7CzrkDUhargrKjimpvLUh3zYd1rCKp6re1inAv88Ejez\nhrMJA5SqHROqBCZ7oC8KcWkMxcKrNgXeXO57wrOOWc8vr6WNbyU8Sy40yMFkYTKEfPuzFFOlBvt9\nUyLoA3XO66ohHZRScKH0rJR3bAuEgs7jr/qll0vHG360AYhK5T/HKcrSw445wBNail5wKhB9s5jH\na1jA2GBXbArcr2BLyWzCgJ+1qERgkl/Tja+xHV99V4u4bTvU9LaxvhubSi+jNtIuegpr4RZLNDqs\nFDgFyhCVaX/2E0U7HjgGGBLZ2NVtBX1LU6EJK6Ugo0HqyZRGtjy4WpswE/0Em/3oVRrRukJjA3GK\nzpqmg92skgGEiScV0KyPEzU84G4aThHFg7ihkhuLSh+hpUCWt3M+oBe1D7MfK2k3FMYgy7ZEFwcR\nbnD4RS5ehXEqsMxHWc7hHzmA9Q36WRT6wU14uAcADC1Qhq2PvHrTKE76aImQN/gG2Fs88GelFGQE\nrwWH/Xp4eF8l6zBhMWPQmxdIgPmhAQlxKeIUZelh5wLPdFGECetinqnhN5qweH+5FhhSwZaSlg9p\n1hz7CC0Fcrif87FgemA/hsOmOGArRTvK11e0YLA04yRreYDdFGg78qZ4HqyIBktGGxon2nDTNnU7\nDgUgCKNkkAwlyzrz4V0NWLcY99ASh+mF98RUWCuFGYnDT913mHBLg8DGK+ANQAIGnChFE8QOBkW9\n20UBanhcr+mJ3DqzmMdruH9X1UM0PVPgbgXbSgoTdpBsPkJLAUw050MBdOLs0DtEmqKjqGdNmgNO\nYeWUlwwWK6h/4bQWA1ww+mIiqjWVWgyH8xV+odlzT4+i0QfIsJLPF7ym6a4lseQlOPcklHNDKdm7\nwLOYnftwmDCnZSXvuE2BFglI0ICji76vKChJnaOlrBYQir+VfL44hp0bvI9p2av6NZpMAG5WDRPN\nwUWy+QQtBeDVnI8Z6+NfXHqHSJNcHv2osQ+ZpsLyF6w5tIyihbV3Bo5q/UDn+Mlg9xXekYfDMNRP\noo0GyNAJ9r7hX2ENQjGYdyulTGMgt98ppl/NIcDrhmWlDT9aTIWlNOq6KPqC3ZP2AKRLYOe5ryg8\nnxp2Ys5kF+TLWgbgZtUwkQ29JJuLaSmAjOZ8LGE6cO3Ve0+WxfSAhzXl31zN5eghN9BY2/eAgaMC\nw1jI/6vGd6UtGJ5Z+15Kg2AZKsHueYQJN03N36DLrBQQ0FTKKIYhNE8Il2340SYgccPirqKkcFrY\nMQ9439RyHzdbSaL8+0g2F9NSoMq2unvVBMPPPr13JRiroPSsfvEObZjZTb94oPd6rosHv3HJI9/V\ngctKG360BUigWxJTzKSwg+7Z3RsnpaWtpHjxjkGyQfQ/RUsB0YrzEUnM8LQJQ+9IgSbTA4Qbr79z\nuJqe3I3bwGBJ5qjcjdCv3/O12bEagFzvUoBWKWHHvOClpKWtpPn6OyTZ4PExWgrIlpyPa5gepC1M\nyUDvWIEm0wPE/MGX0EpwAk7F5NjNEJAuP5IEAsZLaNUK8edoKVBiwfmIJmZsYXWIHGQyPVD433sV\n/BZk9m/tRLvC7Zv56lcQMF4FjyQbOD5ISwHpYrdpNDGDNPN8HCIHmXtMUe4f/EMWD1zr23xyLIWu\nH8t3kkTA+EMWItmASxnO0/HQKzYLKkz4SEon0wPy2iIHhTE9UGEr7rpZgj/448u3pfoPlvXPFUmz\naRRR5DQtZRsjYcK0XwAWRMLZQCDVzfRQetNe/IVAP5VnTUe5+z8et2H77K+4bz8fP4rA3HBnks1h\nWkoYX0MxU0xiRlhKN8Kz3qbASKYHSc5/Ee4GON9NG4H5L8I1UYRM4QAtJWD7LGKhOB8GMSMwpRNJ\nrbchMIDKM+uhpdI7qvRl/pYR+AUEiDoHimqSzVFaStD2WchJcT40MSM0pQtPrTfTAmGVaZfK49jL\nkIdhF8L5XtoIzIOwJtkcp6WEbZ9VnA+TmBGW0gWl1tuijkQyPaRkQfpyZZPvZQQSRcCxA1rzdIJG\nx3hWieJ8xKfcAVELjGV6SMFczdN3Mso/ZwRSQaAGlsvyMHk6caNj4PZZBzMlMOVSVe91LNNDCarp\nXe7qKp8zAqkjwLGpLw5reTSK6Rq8fXbF+QhOudDVexnL9JgFtXl9dMYif/kBBOSbKSxNP01LgcwS\n4HzYTA8DALm9xbiTv2YEkkXA3VyNRf7LR0cJxf2cD5vpYdbQcPt7vUxt8veMwAYC+O848vgHNUug\nDpi72oIAAAAASUVORK5CYII=\n", "text/latex": [ "$$\\frac{1}{2 \\rho_{max}^{2} u_{max} \\left(u_{max} - 2 u_{star}\\right)} \\left(\\rho_{max} \\left(2 u_{max} - 3 u_{star}\\right)^{2} + \\sqrt{- \\rho_{max}^{2} u_{star} \\left(2 u_{max} - 3 u_{star}\\right)^{2} \\left(4 u_{max} - 9 u_{star}\\right)}\\right)$$" ], "text/plain": [ " _____________________________________________\n", " 2 ╱ 2 2 \n", "ρ_max⋅(2⋅u_max - 3⋅u_star) + ╲╱ -ρ_max ⋅u_star⋅(2⋅u_max - 3⋅u_star) ⋅(4⋅u_ma\n", "──────────────────────────────────────────────────────────────────────────────\n", " 2 \n", " 2⋅ρ_max ⋅u_max⋅(u_max - 2⋅u_star) \n", "\n", "______________\n", " \n", "x - 9⋅u_star) \n", "──────────────\n", " \n", " " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_sol[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluating the new flux equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quadratic equations, of course, have two solutions. Here we have to select the positive root, otherwise our model would have an inconsistency. Are you able to see why? To determine its value, we need to fill in some actual numbers into the model.\n", "\n", "Let's start with the same numerical values that we used in [lesson 1](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb) for $\\rho_{\\rm max}$ and $u_{\\rm max}$. But we also have to supply a value for $u^{\\star}$, a quantity that should be experimentally observed in a given road. We propose $u^{\\star} = 0.7\\, u_{\\rm max}$ for this exercise. This would correspond to 84 km/h for a highway with a 120-km/h speed limit, for example.\n", "\n", "Let's numerically evaluate the solutions for $A$ using the following values:\n", "\n", "\\begin{align} \n", "\\rho_{\\rm max} &=10.0 \\nonumber\\\\ u_{\\rm max} &=1.0 \\nonumber\\\\ u^{\\star} &=0.7 \\nonumber\n", "\\end{align}\n", "\n", "Now evaluate the numeric result for each root of $A$ using the `evalf()` function, where you pass the numeric substitution as an argument. It's very cool. Let's try the `[0]`-th solution for $A$, first:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAAAQCAMAAACInqX4AAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAv1JREFUSA2tVtkC4yAIxMQjxmj9/69d\nLkXb7D6tDykIAxMEU4D/urz/r+H+Ecwdx7kmW3SfDXdGkV1oUdzNOiDpirIAfGxthp2ebm6aPX4A\nisMHmITKSx5DBySTLmM99RLz1Sfp8jBp344CsWGaxTohrsvKyLlCaU9F/OKZs4d4U67F3gjD1TEJ\n3vIY2l3EK8+Smu5TiUb66Ey6BfTuuL1YDaJHgBXOhZh1Cm6eRye+14GPxZ5DCw63cHNK8JJnQd/E\nAk4OR9KmG+lYmbTrVLsPnSWepr6SQaSDPujTO7Fu7D48fccTAggPPha7gMhi0kueDc2eqSdC4ZKC\nqj5oQY1iuCnhWMO6Qag/0SFwwzUNK55JSEfaXOxG1aSXPAs6jYwK2PVBC7B7yBHf1seIXc1LrTtk\nrRdc+o5fpCc5tsdy0gTimtJbnkk6wnfGXR+kTzpwTJX6fQJUnq9x6N8h+B5gDmhBb1oa57lJCVIn\nlMSenYczULdP6TWPoXeSFIfLoD+arFJuIc2tKmc7qOwQvzTQPcZb43wI7Z9JWux8cd00nVNK0vJ7\nHkPvGf9CmgKOSpMYO0+7UtlDnDxs5LVcSeoJrn0+Zxz1tyuLQjJjDu6xdHwmex5Dl1FZuXRg1yXZ\nmemL0Vt0OsFRQEplh1yjuuC0N5gIvQStVOHQK0XtchXRdJr0mmdBP5xErhiKuulCqyRaPWPCSyvA\nfEb9Voi+AQZyVIYiMzs8KQE0/jJMO136VOnK179I73nIT9CZjzNSoMLDYDoFYEd+cHVl4+DOntYl\nBHxGw1Y+Ov2Qa5zjoZ6W3pr2xjdHQAImaeg9z4IudBX4C5GFx9N0pPpD2j9YY3/JnTWsK0Q7B2rj\nPyF8MDNOINKBj9bs3BSV7m6TcFp/8yxoSM19Ah23l2GeOrT29LtJs7e7PyTWO0fhvFgNgp9W7hx4\n5E8IH9v0rCEqeLW7HPPN3zaT3vIsaGSb0phc7oMfXXbHs355876FGF/W4b7/voI9Dossk7A8v3l0\n6w/otidHUD5UdgAAAABJRU5ErkJggg==\n", "text/latex": [ "$$0.0146107219255619$$" ], "text/plain": [ "0.0146107219255619" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aval = A_sol[0].evalf(subs={u_star: 0.7, u_max:1.0, rho_max:10.0} )\n", "aval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's the positive root. Let's evaluate the negative root (but we're not going to save this result, because we don't need it!):" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMMAAAAQCAMAAACcCzsZAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAHr8+Xo9uDq9+\n35/P7y5Oa96t3wAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAvpJREFUSA2dVtcW2yAMZTmAJ///tdVC\nkpO08SkPjtYVGggSwn+sWh+CHhs+9PfYrCzLywfp+LqRl7x2XsiJ7EZNiDOsPSV1q5iiQtP3K4RY\n4BOMAqa01DmsL2jQ+9V6CHm1JJSPfVsHWZbBawsmMyoEhZhh7WeI6TgB7yy3rYa+415On9A5Fcuo\nUNMSQ09/Qfv4gS4rCjYu+I2vOXbOQeoBZTWZUc6FGW4RfNWBzs1yGRj+usDH6beWWgERCJUKqQE/\nYPuvaLK2z47G4UXekbrxMwdUhAuLCg3nvBxlEGgpLDQcA5NIw2PqgLpC1w74OD2DUGNUIeCFB2zu\n6NEkts8gYB5ZRDfe4sUjyxYmm9QNIoaNTmcSt2yZOYeOQqe3yI3aMc+5PtFTw795BiD4Oz+jRNu5\ngcmEukPMEDCrhPIWxXTF+h5fOMy4xaSg5rV3mAgWU+e1AopmLUzz0xzo4kDUrxzUEJ2/eB/BHDuy\njbcEivVbqeHVcFKUymMH5ElXgu74jkZftJ7mULW5P3IwQ5iteVMI5sIhqYfmwHq6EnccdKUyjwsf\nuFk1Q8fkVnzchxeNI6b9IwczdJfdxJR0Xa8+u2OXIXqlBMh/haioY33QDH6is1swgPMsycTe+YmG\nkz1L+pnDHWKGRQ7SLe98hkUuK9Hz5YODbpTcQZ2DsygcGss510HB8V2GshuvaIkTDVSmlIeYYcGq\nRB5Lw4As0YsUph4fAfR60nPAVFilD1SGb2jE6NrokHT0G2mujHfxXnqGNXLTOhdBDU9qrPzbkCiW\nA+eBD4jqE91JDQIwSjZZ6JGZO3q0hk9ExOGvKziKNPjGg1orIF1FhMqU8pBpeCb6h0XlVEzDHBo1\n3vR0gk58M4yCwYcO1JWvUdnRoSly98mpXA17VvmaUD6kdIw9cT3nHJrMqGAQePB5Cg7+h0VNVcuz\n9S5xOX3Z+rbTE2tUOPdNTL+hXfhC5jzvBBa88yidz7hA3n8M8m/D820r8lNhVHkZBa/Dp6mI/gCf\nIiXX7Nt+UwAAAABJRU5ErkJggg==\n", "text/latex": [ "$$-0.0171107219255619$$" ], "text/plain": [ "-0.0171107219255619" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_sol[1].evalf(subs={u_star: 0.7, u_max:1.0, rho_max:10.0} )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, numerically evaluate $B$ in the same way:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL8AAAAQCAMAAABwfl4PAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAXo9uDq9+Pt+f\nzx7vLr9ObtN+HgAAAAlwSFlzAAAOxAAADsQBlSsOGwAAA0xJREFUSA29VtuypCAMDMpFVJD//9pN\nJ0DYGetsnZelapxAh6QhFyX6D8P7H538A/5xL0C3bfvqweaLFGpUnfgQXY4f5Pda+0Y3YXJj0cch\npSPqgLOpOY2vMGvskR8yuvTuh3xWrcD66bAD2NykzXkmJjq18cBWHwr/TiwucM6eIhZ9LHTVuzDs\nsIUH7zJNM77ArHzdg3+XXv1cMR+Ntdn4gWfuh1nnhqRNNEUnhxocpnvbiR6QWuCt4TwH6+eLBd9g\nvUeOY2WaZnyBWZe2Nvh36dWPT1dU/mfArl38QrK5SVH4pwp8WKd0cxKlxosGe0yJwk3UGg5QGwdA\ntzyQpiEzvsAwXwb/Ib35gZPOX9VTS1jjYXOTHKhSlGuf/EXbIQgGy3FgOlGQlKzT7oXdpmnGxQ6X\nlPyX2AmQSQA+/GBJ+adhqPOyuUkclBb8oykWrx0l3EeVxQlP/uOUB0dCh64MzdW44H0Dm1NCZBLw\nDz9Y+g1/H7RmeVvmUt6DlHKKp6QLl/GA7xOmQ+fA6cXhkSEdC0Wvmp/8O7xzjulekzhJv/2wzd/w\nL9kdTfuPMD8lj5lOPtFgDH6Q9f4e/E+NGZb0GEPzg3+HC04re02SfV9+ePUX/JE6fmtyt2IvjnLn\nBuNphV19nj32a7eutmugpuYH/w7LrQh/k8QfGtnffgb/S68qNa0fsrlJpG+HdLPOI7mPAtVReeMC\n82Iq3AMRFXIjebijaiCm5mIcmgrvGa+5VqMzCSjGtx+9f7rFslOPrGjzKWnj5NbMRdak68ZWLk31\n2uIKiy+q8k5xuJILXXTeikYBhqbxBb4SRsupmMTvMwn7t5/OP4vNCI8XctvmU7qEDvPnu69y/+Gg\noul0tLTC2438l2AWiah+YDwa5UVzGgf/DkMc9Tuldz8Md/4XXvH+YFp60nU+kCzdreBNJ/lTkD4B\nN8u3xUc2OIB/QEhLlW8erZmo/BdNcyNUevtkefTPKb36kU3QZgrVPQGZ6rVbzLkhFJlM2KT1uBzz\niez3W4ibfqwYXEKMh5C59ZtHM0a+AeDMNM1N/xQBzPE9210lckN69UO13u1URUpJmKkBftr8TfJc\noDr8G1w+jKnqKPd342RwN/3X36sf1fgDHLspZCeSA2oAAAAASUVORK5CYII=\n", "text/latex": [ "$$0.00853892780744381$$" ], "text/plain": [ "0.00853892780744381" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bval = B_sol.evalf(subs={rho_max:10.0, A:aval} )\n", "bval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Turn off $\\LaTeX$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\LaTeX$ output is great when we're looking at algebraic expressions, but it's a little too much for simple numeric output. Let's turn it off for the rest of the exercise:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sympy.init_printing(use_latex=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Green light: take 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's re-examine the green-light problem from [lesson 1](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb) but using our newly computed traffic-flux equation. We shouldn't have to change much—in fact, we can simply create a new `computeF` function and leave the rest of the code unchanged.\n", "\n", "There's one last bit of housekeeping to do so we don't run into trouble -- we used the variables `rho_max` and `u_max` in our SymPy code and we defined them as SymPy `symbols`. You can check on the status of a variable using `type()`. " ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n" ] } ], "source": [ "print(type(rho_max), type(u_max))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we try to use SymPy variables with NumPy arrays, Python is going to be very unhappy. We can re-define them as floats and avoid that messy situation. " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rho_max = 10.\n", "u_max = 1." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def computeF(u_max, rho, aval, bval):\n", " return u_max*rho*(1 - aval*rho-bval*rho**2)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "from matplotlib import pyplot\n", "%matplotlib inline\n", "from matplotlib import rcParams\n", "rcParams['font.family'] = 'serif'\n", "rcParams['font.size'] = 16" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def rho_green_light(nx, rho_light):\n", " \"\"\"Computes \"green light\" initial condition with shock, and linear distribution behind\n", "\n", " Parameters\n", " ----------\n", " nx : int\n", " Number of grid points in x\n", " rho_light : float\n", " Density of cars at stoplight\n", "\n", " Returns\n", " -------\n", " rho_initial: array of floats\n", " Array with initial values of density\n", " \"\"\" \n", " rho_initial = numpy.arange(nx)*2./nx*rho_light # Before stoplight\n", " rho_initial[int((nx-1)/2):] = 0\n", " \n", " return rho_initial" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#Basic initial condition parameters\n", "#defining grid size, time steps\n", "nx = 81\n", "nt = 30\n", "dx = 4.0/(nx-1)\n", "\n", "x = numpy.linspace(0,4,nx)\n", "\n", "rho_light = 5.5\n" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rho_initial = rho_green_light(nx, rho_light)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAECCAYAAAD5OrxGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHxZJREFUeJzt3Xd4VGX6xvHvQ+hRsCNYsVdsu2thdWNDxbYURQTXdf0J\nuoqKBQVcibog2FhREbCuguJiXRXLWqIuroq9IBYUUUERUekh5fn9MUOYiSEkmZNp7/25Lq4Jb055\nrpfhzsk5Z85j7o6IiISjSaYLEBGR9FLwi4gERsEvIhIYBb+ISGAU/CIigVHwi4gEpmm6d2hmun9U\nRKQB3N2i2E5GjvjdPev/DBs2LOM1qE7Vmct15kKNuVRnlHSqR0QkMAp+EZHAKPjXoKioKNMl1Inq\njJbqjE4u1Ai5U2eULOpzR2vdoZmne58iIrnOzPBcvrgrIiKZo+AXEQmMgl9EJDAKfhGRwCj4RUQC\ns9bgN7P2Zva0mVWmoyAREWlctQa/mXUHXgW2AdZ4D6aZrWNmN5vZTDP7yMyeMbNdIq5VREQisLYj\n/kHAYcTCv7b7R6cAnYA93X1X4HWgxMw6RFKliIhEZm3Bf4C7z6ptATM7HDgCuNzdV8SHrwIKgCGp\nlygiIlGqNfjdvS7n9XsAK4H/JqxXBkyLf09ERLJIFHf1dALmunt5tfHZQDsz2yiCfYiISESiCP6N\ngMU1jC+Kv24YwT5ERCQiuo9fRCQwUbReXABsWsN4m/jrj9W/UVxcXPV1UVFRkI9FFRGpTUlJCSUl\nJY2y7To9ltnM7gb+5O6/+g3BzMYBpwGFief5zexxYB9371BteT2WWUSknjL1WOY1pfVDQDOg86oB\nM2se//tDDS9NREQaQ32Cv8afNO7+H+AZ4CozaxUfHgqUASNSK09ERKJW6zl+M7sGOBzYEnAze4fY\nkf++8Xv1V+kJjALeNbMK4GugyN3nNU7ZIiLSUGq9KCKSA9R6UUREGkzBLyISGAW/iEhgFPwiIoFR\n8IuIBEbBLyISGAW/iEhgFPwiIoFR8IuIBEbBLyISGAW/iEhgFPwiIoFR8IuIBEbBLyISGAW/iEhg\nFPwiIoFR8IuIBEbBLyISGAW/iEhgFPwiIoFR8IuIBEbBLyISGAW/iEhgFPwiIoFR8IuIBEbBLyIS\nmEiC38x+Y2ZPmdkMM3vfzF43s55RbFtERKKVcvCb2dbA88B8YDd37wTcCfzLzI5JdfsiIhKtKI74\nuwLrAje4eyWAu48HFgEnR7B9ERGJUBTBXx5/bbZqwMwMKIho+yIiEqEogvl+YCZwmZkVmlkTYAix\nHwTjIti+iIhEqGmqG3D3xWZ2KHA3sABYAvwCHO7ur6S6fRERiVYUF3d3BKYDXwLru/vGwFDgETM7\nMtXti4hItFI+4geuAtoA57l7KYC7P2BmvYB/mlkHd69IXKG4uLjq66KiIoqKiiIoQ0Qkf5SUlFBS\nUtIo2zZ3T20DZh8D7u67VBsfBVwMbO/usxLGPdV9ioiExsxwd4tiW1Fc3P0e6GBmBdXGtwIqgZ8i\n2IeIiEQkiuC/idipnitXDZjZwUA34AF3XxjBPkREJCIpn+oBMLMuwKXApkAFsSP9e4Ax7l5WbVmd\n6hERqacoT/VEEvz12qGCX0Sk3rLtHL+IiOQQBb+ISGAU/CIigVHwi4gERsEvIhIYBb+ISGAU/CIi\ngVHwi4gERsEvIhIYBb+ISGAU/CIigVHwi4gERsEvIhIYBb+ISGAU/CIigVHwi4gERsEvIhIYBb+I\nSGAU/CIigVHwi4gERsEvIhIYBb+ISGAU/CIigVHwi4gERsEvIhIYBb+ISGAiC34z62FmL5nZm2Y2\ny8ymm1nfqLYvIiLRiCT4zWwgMATo7e6/AXYEPgUOiWL7IiISHXP31DZgtjUwE+js7m8ljLcHOiSO\nxcc91X2KiITGzHB3i2JbTSPYxinAT9UD3t3nAfMi2L6IiEQoilM9BwBfxc/xv2xmH5vZNDM7LYJt\ni4hIxKI44t8C2Aq4AOjm7vPNrCdwv5m1d/cREexDREQiEsU5/s+BbYDfu/urCeMPAkcCG7v78oRx\nneMXEamnbDvHvxhw4N1q4+8C3YGdgbcTv1FcXFz1dVFREUVFRRGUISKSP0pKSigpKWmUbUdxxH8/\n0Ato4+5LEsYHA8OB37n7mwnjOuIXEamnKI/4o7i4++/46x7VxncDlgEfRbAPERGJSBTB/wAwHfi7\nmRUCmNmBQA9geOL5fRERybyUT/UAmNn6wCigC7ACKAXGuPsdNSyrUz0iIvUU5ameSIK/XjtU8IuI\n1Fu2neMXEZEcouAXEQmMgl9EJDAKfhGRwCj4RUQCo+CXILk7387/icrKykyXIpJ2Cn4Jzoeff0vn\n069h866X0uWcGzNdjkjaRfGQNpGcsHzFSq6640muvedZyitiR/rPvzGTuT/8TIeN18twdSLpo+CX\nIDz/xsf0HzGJWd/88KvvLV66AjbOQFEiGaLgl7z2w0+LuXD0g9w79bU1LrNsxco0ViSSeTrHL3nJ\n3fnnE/9j557DkkK/7TqtGD+kD/vt3rFqbKmCXwKjI37JO5/N+Z4zr57EC9M/SRo/8fB9+MeFvWi/\nUVumPLe6N5CO+CU0Cn7JGyvLyrnmn8/w9zunUrqyvGp8q/YbMvaS3nT9/e5VY4Wtmld9reCX0Cj4\nJS9Me/dz+o2YyIwv5lWNNWlinN/7UK488zgKW7VIWr51y9XBv3R5adrqFMkGCn7JaT8vXsalNz3C\n+IdfThrfZ+ctuW3oKey105Y1rpf4g0BH/BIaBb/kJHdnynNvcd51D/Ddj4uqxgtbtWD4X4/nnBMP\npqBgzfcu6IhfQqbgl5zz1bwfOXvU/Tz53w+Sxo89sBM3X9KbLTfdYK3b0Dl+CZmCX3JGeXkFYya/\nwN/G/TsprNtv1JabLj6J7ofshVndGhS1bpF4xK/gl7Ao+CUnvPXxV/QbPpG3Z86pGjMzzup5ECPO\n7kbbdVrVa3uJp3p0xC+hUfBLVluybAWXj/s3N05+gcrK1b2ad9u2AxOG9mX/Tts2aLu6uCshU/BL\n1nrilfc5e9T9zPluYdVYyxbNuPz/juaiU7rQrGlBg7eddHF3hS7uSlgU/JJ15i34hfOue4Apz72V\nNH7o73Zi3OA+bLfFJinvQxd3JWQKfskalZWVTHj4FS69+RF+WbK8anzDtoWMvuBE+nbdt84Xb9cm\n+XZOBb+ERcEvWeGjWXPpN3wir74/K2n81GP257rze7LReutEuj9d3JWQKfglo1aUlvH3O57kmnue\npay8omp8uy02Ydzgkzn0dzs3yn4TL+7qA1wSmsiD38xeAToDW7v7nLUtL+F6YfpM+o+YxOdfz68a\na9a0gEtOPYKhf+lKyxbNGm3fOuKXkEUa/GbWg1jo+9qWlXAt+HkJF46ewj1PJjdH6bzHtowf0pdd\nt+3Q6DXodk4JWWTBb2bNgZHAVKBrVNuV/OHu3Pvka1wwego//rK0arztOq245tzu/N8ff0+TJunp\nDZR8O6eCX8IS5RH/2cDrwGco+KWaz+Z8z1kj7+P5N2YmjSc2R0mnQp3qkYBFEvxmtgFwEbAfcHoU\n25T8sLKsnOvufZYrb38yqTnKlptuwNhLT+bohOYo6ZR4/WBFaRkVFZW1Ps1TJJ9EdcR/OXCvu38d\n1X3WkvtefW8W/YZP5KMv5laNrWqOckX/Y1mndcuM1WZmtG7ZvOpof3npyozWI5JOKQe/mW0PnADs\nlHo5kg9+XryMwTc/wriHkpuj7L3TlkwY2pd9dt4qQ5UlK2zVoir4l61Q8Es4ojjiHwVc7e6L67pC\ncXFx1ddFRUUUFRVFUIZkmrvz0PNvc+51DzBvwS9V44WtWnDVmccxoNfBNE3h+TpR06d3JZuVlJRQ\nUlLSKNs294bfeWlmBwK3A7u6e3l8rJjYqZ+O7v5VDet4KvuU7DTnu4WcPeo+nngluTnKMQfuzs2D\nerNV+w0zVNma7XpicVWP3g8fGJaW20hFGsrMcPdIzqWnesR/GFAATE84t79p/HWqma0EBrv70ynu\nR7JURUUlNz3wApfd+u+kT8BuumEbbrr4JHocundkz9eJmtovSqhSCn53HwYMSxwzs1VjR+mTu/nt\n7Zlz6Df8Xt76OPmf+cweB3H1Od1Yb93WGaqsbvTpXQlVYzyrx6q9Sp5ZsmwFw8Y/zj/ufz6pOcqu\n28SaoxywR8Oao6RbYcuE5/Uo+CUgUX5y9yhgBLFTPU7sVE+pu+8d1T4k85787wf8deR9Sc1RWjRv\nyt9OP5qL/9SF5s1y57l/OuKXUEX2v9TdnwKeimp7kl3W1BzlkN/uyLjBfdh+y3YZqqzh1IxFQpU7\nh2eSEZWVldz2yH+55KaHf9Uc5YaBJ3DK0ftl7cXbtdHFXQmVgl/WaMYXseYo095Lbo7yp6P34/qB\nJ0TeHCXd9IROCZWCX35lRWkZw++cyqh/PpPW5ijppiN+CZWCX5K8+OYn9B8xkc/mrG6O0rSgSVVz\nlFYJYZnrdHFXQqXgFwB+/HkJF934IHc//r+k8f07bcOEIX3ZbbvNMlRZ49GjmSVUCv7AuTsTp77O\nBaOnsODnJVXjbQpbMmpAd/p1PzBtzVHSTc/qkVAp+AM265sfOHPEJJ574+Ok8Z6H7s2NF/Wiw8br\nZaiy9NDFXQmVgj9AZeUVVc1RVpSWVY1v0W59xl56Mscc2CmD1aVPcvtFXdyVcCj4A/O/92PNUT6c\nldwc5byTDuHKM48L6pn0urgroVLwB+KXJcsZcssj3PrgyyQ+FjvbmqOkU+Ind3WOX0Ki4M9z7s7D\nL7zDgGsnJzVHad2yOX8/6/isa46STjril1Ap+PPYnO8Wcs6o+3n8lfeTxrt23o2xl56clc1R0kkX\ndyVUCv48tKbmKO02bMOYi3pxwmH75OzzdaKkT+5KqBT8eeadmXPoN2Iib85I7nrZv/tBjByQ/c1R\n0klH/BIqBX+eWLq8lOIJjzP6vuepqKisGt9lm/aMH9KX3++5XQary07Jt3OuxN31m5AEQcGfB56a\n9iFnjbyPr+b9WDXWonlTLvtLVwadekRONUdJp2ZNC2ha0ITyikoqKiopK6/QXEkQ9C7PYd//uIjz\nr/8Xk5+dnjR+8G9izVF22Cr3mqOkW2GrFlV9BpYuL1XwSxD0Ls9BlZWV3PHYNAaNeZifFy+rGt+g\nbSHXn9+TU4/ZX6cs6qh1y+ZVwb9sxUrWb1OY4YpEGp+CP8d8/OU8+g2fyH/f/Txp/JSu+3H9wJ5s\nvP66GaosN+kCr4RIwZ8jVpSWMeKupxh599NJzVG23Xxjxg3uw2H75kdzlHTTEzolRAr+HFDy5if0\nHzGJT+d8XzXWtKAJg/50BJednl/NUdJNDdclRAr+LPbjz0u4+MaHuOvxV5PG87k5SrrpQ1wSIgV/\nFnJ3Jj31OgNv+HVzlJEDutM/j5ujpFvrFjril/Ao+LPMrG9+4KyrJ/Gf15Obo/Q4ZG/GXJz/zVHS\nTRd3JUQK/ixRVl7B9RP/wxW3PfGr5ii3XNKbYw/aI4PV5a/qn94VCUEkwW9mewJnA52BcqAAeA64\nyt0XRLGPfPbaB1/Qb/hEPvj826qxJk2MAb0O5qozj2fdwnCao6SbLu5KiKI64p8MfADs4+7LzawD\n8DxwpJnt4e4rItpPXlm0ZDlDbnmUsQ++lNQcZc8dtmDC0L78dtetM1dcIHRxV0IUVfBXApe4+3IA\nd59rZtcCtwNdgYcj2k/eeOTFdzjnmsnM/eHnqrHWLZtz5ZnHcd5JhwTbHCXd1IxFQhRV8Hdy9/Jq\nY/Pir7oameDr7xYy4NrJPPbSe0njRx2wG2Mv7c3WHTbKUGVhSry4qyN+CUUkwV9D6APsADjwchT7\nyHUVFZXcMqWEoWMfZcmy5OYoN154Iice/hs9XycDko/4y2pZUiR/NMpdPWZWAJwO3O7un69t+Xz3\n7idf02/4RKbPmJ003q/bgYwc0E0PBsugQp3qkQA11u2cfwNKgfMbafs5YenyUq6Y8AQ33PdcUnOU\nnTu2Z8JQNUfJBrq4KyGKPPjN7DSgJ1C06mJvdcXFxVVfFxUVUVRUFHUZGff0q7HmKLPnrm6O0rxZ\nUy47vSuD/tSFFs2bZbA6WSXpA1ylOuKX7FFSUkJJSUmjbNsSbyNMeWNmpwCDgEPdff4alvEo95lt\n1tQcpWifHRg/pK+ao2SZl976lKL+1wNw4F7b8fJtF2e4IpGamRnuHsmFwMiO+M2sL9VC38yOAdq7\n+21R7SdbVVZWcue/X+XiGx/6VXOU687rwZ+PPUAXb7OQbueUEEX1yd0+wG3AZUCXhIA7EJgbxT6y\n2cdfzqP/iIm88k7ydew+R/2O0RecqOYoWSz5k7u6q0fCENUR/xigOXBttXEHrohoH1mndGUZV9/1\nNCPueiqpOco2m23EuMF9OHy/XTJYndSFLu5KiKK6j3/DKLaTS15661P6j5jIJ18lN0e5sO/hXH7G\nMUmBItlLT+eUEOnpnPW08JelDBrzEHc8Ni1pfL/dOzJh6CnsruYoOUVH/BIiBX8duTv3PzOdgTf8\ni/kLF1eNtylsydXndKN/94MoKFBzlFzTqsXq22qXl5ZRWVmpJjeS9xT8dfDFNz9w1sj7ePa1GUnj\n3Q/eizEX92KzTdbPUGWSqiZNmtCqRTOWx3sgLC8tSzr9I5KPFPy1KCuv4IZ4c5TlCc1RNm+3PrcM\n6s1xf1BzlHzQumXzqn/fZStWKvgl7yn41+D1D7+k3/CJvP/ZN1VjZsa5J6k5Sr4pbNWCH39ZCugC\nr4RBwV/NoiXLGTr2UW6ZktwcZY8dNue2oaeoOUoe0gVeCY2CP8GjJe9yzjX38+385OYoV/Q/lvN7\nH6rmKHlK7RclNAp+4Jvvf2LAtZN5tOTdpPEjD9iVsZecTMfN1BwlnyUf8Sv4Jf8FHfwVFZWMnVLC\n0FsfY/HS1W2BN9lgXW68sBe9uqg5Sgj0vB4JTbDB/96nseYob3w0O2n8jG6/Z9SA7mqOEpDClmq/\nKGEJLviXrVjJFRMe5/pJyc1Rdtp6UyYM7cuBe22fweokE3TEL6EJKvif+d9HnDXyPr78dkHVWPNm\nTRn6l6O45NQj1BwlUHpej4QmiOCfv3ARA2+Ywn1Pv5E0/oe9d2D8kD7suPWmGapMskHrlqt/4C9V\n8EsA8jr43Z07H5vGxWMe4qdFq5ujrN+mNded15PTjlNzFNGpHglP3gb/zNnf0X/ERF5++7Ok8T5H\n/Y4bBp7AJhu0yVBlkm0ST/Xo4q6EIO+Cv3RlGSPvfpoRdz3NyrLyqvGO8eYoXdQcRarREb+EJq+C\n/+W3P6X/iEnMnP1d1VhBQRMuUnMUqYUu7kpo8iL419QcZd/dOjJhaF86bb95hiqTXKBP7kpocjr4\n3Z3Jz0zn/GrNUdYtbMnws47nrycUqTmKrFWhTvVIYHI2+L/8dgFnjZzEM/9Lbo7S7eA9GXPRSWze\nTs1RpG6SjvhX6OKu5L+cC/6y8gr+cd9zDBv/eFJzlM02WY+bB/Xmj0V7ZrA6yUW6uCuhyangn/7R\nbM4Yfi/vfZrcHOXsE/7A8L/+kTbrtMpgdZKrkm/nVPBL/suJ4F+8dAVDxz7Kzf8qSWqO0mn7zZkw\ntC/77tYxg9VJrtMRv4Qm64P/sZJ3OefayXzz/U9VY61aNKO437EM7HMYzdQcRVKkRiwSmqwN/m/n\nx5qjPPJicnOUI/bfhVsv7aPmKBIZtV6U0GRd8FdUVHLrgy8xZOyjv2qO8o8LTuSkI36r5+tIpPQB\nLglNJMFvZpsAo4F94kMfAOe7+7f12c77n31Dv+ETef3DL5PGTz++M9ec24MN2qo5ikSvWdMCCgqa\nUFFRSVl5BWXlFTqFKHkt5eA3s+bAf4CZwKoH4dwJvGhme7n70rVtY9mKlVx52xNcP/E/lCc0R9lx\nq3ZMGNqXg/beIdUyRdbIzGjdsnnVb5jLVqykre4QkzwWxRH/qcDuwPHuXglgZpcA3wJnAdfVtvKz\nr83gzKsn/ao5yuA/H8ng045UcxRJi0IFvwQkiuDvAXzl7rNXDbj792Y2I/69GoN//sJFXDB6CpOe\nSm6OctDe2zNucB927tg+gtJE6kYXeCUkUQR/J2KneaqbDRxS0wqrmqMs/GX1WaD127Tm2nN7cNpx\nB9CkiZ6vI+mlC7wSkiiCfyNgcQ3ji4DWZtbC3ZMOoU6/6p6kBXsf8VtGX3Ai7TZUcxTJjMQj/klP\nvc6r783KYDUijSuK4Pe1L1KzjpttxNhLenPkAbtFUIZIwyUG/zX3PJvBSkQaXxTBvwBYt4bxNsDS\n6kf7ADbvLfbvtA1Fu29Cy5ULalhVJL323GFzXnzzk0yXIbLa4rmxP43AEp9906ANmD0F7OTuHauN\nfwAsdvcDqo37OzPnsOeOW6S0X5Eo/bRoKeMffoXZc3/MdCkiNRo/tC/uHsmnV6MI/jOA8UBHd/8q\nPtYO+Aa41N2vr7a8p7pPEZHQmFlWBX8z4E3gY6APsXP+dwAHAHu5+7Jqyyv4RUTqKcrgT/m+SXcv\nAw4HKoAZ8T/rAIdUD30REcm8lI/4671DHfGLiNRbVh3x56uSkpJMl1AnqjNaqjM6uVAj5E6dUVLw\nr0GuvBlUZ7RUZ3RyoUbInTqjpOAXEQmMgl9EJDAZubib1h2KiOSJrLmPX0REcotO9YiIBEbBLyIS\nGAV/ljKz9mb2tJlVrn3pzMmVOkVygZm9YmaVZrZlY+4niscyA2BmmwCjgX3iQx8A57v7t3VYtxlw\nOdATKCfWxGWQu0+Lqr6I6pwN/FTDty509xcirLE7cD1QRj37HaR5LlOpczbpmcs9gbOBzsTmowB4\nDrjK3Wt9Jnia5zKVOmeThrmM72tbYr20i+JD6wLfAyPdfepa1k3nfKZS52zSNJ/V9tuD2L9/nf4v\npTSf7p7yH6A58B7wALHfIpoAdwOfAoV1WH8csfaNG8b/fjqwFNgjivoirPPLKOupZT+vAdvGa6us\n57ppmcsI6kzXXM4EpgCt4n/vQOyBgp8ALbNoLlOpMy1zGd/XOcDXwDbxvxswIh48B2XRfKZSZ9rm\nM2GfzYHPgCeASmDLOqzT4PmMqugz4sVunTDWLj7JF61l3R2JPeDtz9XGPwSeiHhyG1xnOt8QQJP4\na70CNZ1zmUqdaZ7LGav+8yeM/SX+PuieRXPZoDrTOZfxff0R+Eu1sbbxOq/LovlsUJ3pns+EfQ4E\nJgLD6hL8qc5nVOf4ewBfufvsVQPu/n38zdxjLet2I/bT+MVq4y8CXcysdUQ1plpn2rh7Q8+Xp3Mu\nU6kznTq5+xfVxubFX9erZb20ziUNrzOt3P1Rd7+z2nDb+OsPtaya7vdmQ+tMOzPbALgIGExsjuoi\npfmMKvg7AV/WMD4b2L0O61YAc6qNf0nsGsQuqRZXbV8NrRMAMxtlZtPM7BMze8bMjo2wvlSlcy5T\nlo65dPfyGoZ3IHYe9eVaVk3rXKZQJ5C596WZbQbcArwVf12TjL4361HnquXTOZ+XA/e6+9f1WCel\n+Ywq+DcCFtcwvghobWYt1rLuMo//nlJtXYANI6gvcV8NrRNgPvCWu3cGdgUeAx4zs7MjrDEV6ZzL\nVGVkLs2sgNi50Nvd/fNaFs3oXNajTsjAXJrZtmb2ObHz6AZ0c/cltaySkflsQJ2Qxvk0s+2BE4Dh\n9Vw1pfmMKvhz5eO/KdXp7vu6+7/iX5e7+1hgKjCiDj80JEEG5/JvQClwfiPuIwp1rjMTc+nus9x9\nO2KnTz4D3jOzzo2xr1Q0pM40z+co4Gp3r+mAtNFEFfwLiN0uVV0bYKm7l65l3UIzq35uq038Ncru\n16nUuSZvxLeZDadR0jmXjaFR59LMTiN269tR7r58LYtnbC7rWeeapOV96e6L3X0gsVslx9ayaEbf\nm/Woc00in08zO5DYbxTjavr2WlZPaT6jCv73gY41jHckdp98bd6L17FFDeuWEbvwGpUG12lmLc2s\nsIZvVcRfC1KsLQrpnMsGy8RcmtkpwAXEWoLWel98XEbmsr51pnsuzazlGr71IbBb/N7ymqR1Phta\nZ5rn87D49qab2Ttm9g7QP/69qfGxI9ewbmrzGdGtSKtuk9wqYaxdvIALqy3bjvjD4eJ/34HYpJ5a\nw21Jj0dRX0R1/hkYV8M2HwaWAa2jrDW+7buBilq+n7G5TLHOtM4l0JfYD/ZNEsaOAc7IprlsYJ3p\nnssSYL8axt8AFmbLfKZQZ9r/n1fbT423c0Y9n1EV24zYT6DJxH6CNQHuIvbhk9YJy3WOFzu22vq3\nkvxBhNOIfRChU8ST2uA642+IX4DfJIz1ii9X3EhvgrtZw/3xmZ7LVOpM51wCfYDlwIXEgnXVn/HA\nsGyZy4bWme73JbHbBZ8BNoj/3YBziYVVcRbNZ4PqzMT/82p1F/Prg9PI5zOSRza4e5mZHU7sUQgz\niF1E/YDYr6vLEhZdTOyj0HOrbWIAsZ9008ysjNiV6S7u/n4U9UVU51Rgc2Bs/NfE9YCFQH93vz3K\nOs3sGuBwYEvA478COrCvu5fVUiOkaS5TrDNtcwmMIfapyGurjTtwRfzrJTXUCGmcyxTqTOdcAgwF\n/g94yczKgZbEzjf3cff748tk/L2ZQp3pnk8AzOwoYp8s3pTYv/lUMyt1971phPennscvIhIYPZ1T\nRCQwCn4RkcAo+EVEAqPgFxEJjIJfRCQwCn4RkcAo+EVEAqPgFxEJjIJfRCQwCn4RkcD8P83i6WRx\nWAVCAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pyplot.plot(x, rho_initial, color='#003366', ls='-', lw=3)\n", "pyplot.ylim(-0.5,11.);" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def ftbs(rho, nt, dt, dx, rho_max, u_max):\n", " \"\"\" Computes the solution with forward in time, backward in space\n", " \n", " Parameters\n", " ----------\n", " rho : array of floats\n", " Density at current time-step\n", " nt : int\n", " Number of time steps\n", " dt : float\n", " Time-step size\n", " dx : float\n", " Mesh spacing\n", " rho_max: float\n", " Maximum allowed car density\n", " u_max : float\n", " Speed limit\n", " \n", " Returns\n", " -------\n", " rho_n : array of floats\n", " Density after nt time steps at every point x\n", " \"\"\"\n", " \n", " #initialize our results array with dimensions nt by nx\n", " rho_n = numpy.zeros((nt,len(rho))) \n", " \n", " #copy the initial u array into each row of our new array\n", " rho_n[0,:] = rho.copy() \n", " \n", " for t in range(1,nt):\n", " F = computeF(u_max, rho, aval, bval)\n", " rho_n[t,1:] = rho[1:] - dt/dx*(F[1:]-F[:-1])\n", " rho_n[t,0] = rho[0]\n", " rho_n[t,-1] = rho[-1]\n", " rho = rho_n[t].copy()\n", "\n", " return rho_n" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sigma = 1.\n", "dt = sigma*dx/u_max\n", "\n", "rho_n = ftbs(rho_initial, nt, dt, dx, rho_max, u_max)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from matplotlib import animation\n", "from JSAnimation.IPython_display import display_animation" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig = pyplot.figure();\n", "ax = pyplot.axes(xlim=(0,4),ylim=(-1,8),xlabel=('Distance'),ylabel=('Traffic density'));\n", "line, = ax.plot([],[],color='#003366', lw=2);\n", "\n", "def animate(data):\n", " x = numpy.linspace(0,4,nx)\n", " y = data\n", " line.set_data(x,y)\n", " return line,\n", "\n", "anim = animation.FuncAnimation(fig, animate, frames=rho_n, interval=50)\n", "display_animation(anim, default_mode='once')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That definitely looks different! Do you think this is more or less accurate than our previous model?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Dig Deeper" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The new traffic-flux model that we developed here changes the way that traffic patterns evolve. In this lesson, we only experimented with the most basic scheme: forward-time, backward-space. Try to implement the green-light problem using one of the second-order schemes from [Lesson 2](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_02_convectionSchemes.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Neville D. Fowkes and John J. Mahony, *\"An Introduction to Mathematical Modelling,\"* Wiley & Sons, 1994. Chapter 14: Traffic Flow." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "###### The cell below loads the style of the notebook." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "css_file = '../../styles/numericalmoocstyle.css'\n", "HTML(open(css_file, \"r\").read())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.3" } }, "nbformat": 4, "nbformat_minor": 0 }