{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solutions Q 8 - 12" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# import all python add-ons etc that will be needed later on\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sympy import *\n", "from scipy.integrate import quad,odeint\n", "from scipy.optimize import fsolve\n", "init_printing() # allows printing of SymPy results in typeset maths format\n", "plt.rcParams.update({'font.size': 16}) # set font size for plots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q8 answer\n", "(a) The variables in the equation can be separated and then integrated giving, $\\displaystyle \\ln(y)=\\int \\sin(t^2)dt+c$. \n", "\n", "However, the right-hand side cannot be integrated analytically but is known in terms of a Fresnel function, which is itself an integral and which can only be evaluated numerically. You could now integrate this integral using Simpson's rule, for example, over the range $-8 \\to 8$ then use $t_0 = 0$, $y(t_0) = 3$ used to determine the integration constant $c$ which is $\\ln(3)$.\n", "\n", "(b) The results of the calculation using the Euler method with the code below, together with the solution from part (c), is shown in Fig. 31. The subroutine is modified slightly to allow the time increment to be negative to allow calculations when $t \\lt 0$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#-----------------------------\n", "def Eulerf(f,t0,y0,maxt,N,s):\n", " \n", " Eulery= np.zeros(N,dtype=float)\n", " dtime = np.zeros(N,dtype=float)\n", " h = (maxt-t0)/N\n", " if s < 0:\n", " h = -h\n", " y = y0\n", " t = t0\n", " Eulery[0]= y0\n", " dtime[0] = t0\n", " for i in range(1,N):\n", " y = y + h*f(t,y)\n", " t = t + h\n", " Eulery[i] = y\n", " dtime[i] = t\n", " pass\n", " \n", " return Eulery,dtime\n", " #---------\n", "\n", "# sample calculation remove hashes to see plot\n", "\n", "dydt = lambda t,y : y*np.sin(t*t)\n", "\n", "t0 = 0.0\n", "y0 = 3.0\n", "maxt= 7.0\n", "N = 500\n", "\n", "soln0,time0 = Eulerf(dydt,t0,y0,maxt,N,0) # positive t\n", "#plt.plot(time0,soln0,color='red',linewidth=3)\n", "\n", "soln1,time1 = Eulerf(dydt,t0,y0,maxt,N,-1) # negative t \n", "#plt.plot(time1,soln1,color='red',linewidth=3)\n", "\n", "#### now do same calculation using built in differential equation integrator odeint\n", "#-------------------------\n", "def dY_dt(Y, t): # returns dY/dt, etc only used in odeint()\n", " dAdt = Y*np.sin(t**2) # equation here \n", " return dAdt\n", "#--------------------------\n", "C0 = 3 \n", "numt = 500\n", "tp = np.linspace( 0,10, numt ) # time as initial, final then number of points\n", "\n", "solnp = odeint(dY_dt, C0, tp ) # solve equations numerical integration \n", "\n", "#plt.plot(tp,solnp,color='grey',linewidth=1)\n", "\n", "tn = np.linspace( 0,-10, numt ) # now do same for negative time\n", "solnn = odeint(dY_dt, C0, tn ) # solve equations \n", "#plt.plot(tn,solnn,color='grey',linewidth=1)\n", "#plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![Drawing](num-methods-fig31.png)\n", "\n", "Figure 31. The Euler method with $N = 500$ over the range $-8 \\to 8$ (red line). The faint grey line is the numerical solution from part (c) also with initial value $y(0) = 3$.\n", "_____\n", "## Q9 answer\n", "Use Algorithm 13 with the potentials $U(r)=\\alpha/r$ or $U(r) = \\alpha e^{-\\beta r}/r$, and initially choose $\\alpha$ and $\\beta$ to be unity. Do not forget to change the force as necessary in the calculation. In other respects, the calculation is similar to that in the example. As these potentials are repulsive when $\\alpha >0$ the trajectories should be scattered like a hard sphere. When $\\alpha \\lt 0 $ the Yukawa the potential is both attractive and repulsive, but in a gentler way than for the L-J potential. The result is shown in the next figure when $\\alpha=-1,\\; \\beta = 1$. The grey-scale shows the magnitude of the potential as a 'density'. The circle represents where the potential has its minima. \n", "\n", "![Drawing](num-methods-fig31a.png)\n", "\n", "Figure 31a. Scattering from a Yukawa potential $U(r) = \\alpha e^{-\\beta r}/r$ with $\\alpha = -1,\\, \\beta = 1$. The impact parameter varies from $d = 0$ and increments in steps of $0.25$.\n", "____\n", "\n", "## Q10 answer\n", "Using the result of question 7, $\\displaystyle \\chi=\\pi-2\\tan^{-1}\\left( \\frac{1}{\\sqrt{d^2/b^2-1}}\\right)$\n", "\n", "differentiation wrt. $b$ can be performed directly, using the standard result \n", "\n", "$$\\displaystyle \\frac{d\\tan^{-1}(x)}{dx}=\\frac{1}{1+x^2}$$\n", "\n", "or by using Sympy gives, after considerable algebraic simplification \n", "\n", "$$\\displaystyle \\frac{d\\chi}{db}=-\\frac{2}{\\sqrt{d^2-b^2} }$$ \n", "\n", "The differential cross section becomes $\\displaystyle I(E_0,\\chi) = \\frac{b\\sqrt{d^2-b^2}}{2\\sin(\\chi)}$. \n", "\n", "In the next step substitute for $b$ using the equation for $\\chi$ above. This produces a difficult looking equation \n", "\n", "$$\\displaystyle I(E_0,\\chi)= \\frac{d^2\\sin(\\chi/2)\\cos(\\chi/2)}{2\\sin(\\chi)}$$\n", "\n", "The trig terms simplify considerably (using the standard form $\\sin(2x)=2\\sin(x)\\cos(x)$) and this produces the nice result \n", "\n", "$$\\displaystyle I(E_0,\\chi)=\\frac{d^2}{4}$$\n", "\n", "As expected by intuition, the result does not depend on the angle $\\chi$ because the interaction between the particles is that of a hard sphere.\n", " \n", "## Q11 answer\n", "Use the Verlet algorithm with values as in equation 30. The initial velocity is +10 and is positive so that the ball is thrown upwards." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Question 11 Verlet method, ball bouncing under gravity\n", "\n", "y0 = 30.0 # metres initial conditions\n", "v0 = 10.0\n", "c = 0.5 # damping constant, 1/ second; see question\n", "t = 0.0\n", "dt = 0.005 # seconds\n", "g = 9.8 # m/s/s\n", "n = 2500\n", "height= np.zeros(n,dtype=float)\n", "atime = np.zeros(n,dtype=float)\n", "velo = np.zeros(n,dtype=float)\n", "height[0]= y0\n", "atime[0] = t\n", "velo[0] = v0\n", "accln = -(g+c*v0)\n", "v1= v0 + accln*dt\n", "y = y0\n", "yold = y0 - v0*dt\n", "for i in range(1,n):\n", " ynew = 2*y - yold + accln*dt**2\n", " v = (ynew - yold)/(2*dt)\n", " yold = y\n", " y = ynew\n", " if y < 0: # bounce ?\n", " y = -y # reverse values on bounce\n", " yold = -yold\n", " v = -v\n", " pass\n", " height[i]= ynew # save values\n", " velo[i] = v\n", " atime[i] = t\n", " t = t + dt\n", " accln = -(g+c*v)\n", " pass\n", "\n", "#plt.plot(atime,height,color='blue')\n", "#plt.plot(atime,velo,color='red')\n", "#plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![Drawing](num-methods-fig32.png)\n", "\n", "Figure 32. Height and velocity of a ball initially thrown upwards and then bouncing but which suffering air resistance proportional to its velocity. The equation of motion is $\\displaystyle m\\frac{d^2y}{dt^2}+mc\\frac{dy}{dt}+mg=0$. $c$ is the damping coefficient with units of s${-1}$. The equation is solved using the Verlet method. If the initial height is large enough or the drag coefficient $c$ is increased then the ball's velocity becomes constant before bouncing due to the effect of the air resistance.\n", "______\n", "\n", "The height profile is what experience would dictate. The velocity initially decreases, reaching zero when the ball is at its highest point and then increases in a negative sense. However, the viscosity of the air offers resistance to the ball and its velocity reaches a maximum downwards (negative) value just before the ball hits the ground. The acceleration of the ball, which is the gradient of the velocity, is always decreasing. After bouncing, the acceleration is always changing and never reaches a constant value.\n", "\n", "**Exercise:** Consider the problem of a stone being thrown upwards at $10$ m/s from a height of $50$ m and falling vertically into water. In air, the friction coefficient is $1$, but $100$ in the water. Calculate what happens.\n", "\n", "## Q12 answer\n", "\n", "The answer is up to you to check!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.6" } }, "nbformat": 4, "nbformat_minor": 2 }