{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary Value Problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we set up the boundary conditions for our differential equations.  We have $\\vec x(0)=[0,0]^T$ and $\\vec x(5)=[j,h]^T=[100,10]^T$.  We also set the step size for our approximations with $h=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set the values for the boundary conditions\n",
    "j=100.0\n",
    "h=10.0\n",
    "T=5.0\n",
    "g=-9.8\n",
    "step=0.1 # This is the step size used in the approximation h=0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the analytic solution to the differential equations when the acceleration is constant.  The acceleration is $0$ in the $x$ direction and is $g$ in the $y$ direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define a function to compute the x and y positions at time t with initial velocity vx or vy\n",
    "def x(vx,t):\n",
    "    return vx*t\n",
    "def y(vy,t):\n",
    "    return 0.5*g*t*t+vy*t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compute the appropriate vx and vy to satisfy the boundary conditions\n",
    "# BEGIN STUDENT: supply the initial velocities in the x and y directions\n",
    "vx=\n",
    "vy=\n",
    "# END STUDENT\n",
    "\n",
    "# This uses the analytic solution with the initial velocities to compute the trajectory of the projectile\n",
    "t=np.arange(0,T+step,step)\n",
    "X=x(vx,t)\n",
    "Y=y(vy,t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we plot the analytic solution using the initial velocity supplied above.  The start and end of the trajectory are marked by red circles.  The solution is correct if the blue line starts and ends on the circles, which indicates that the projectile starts at $(0,0)$ and ends at $(100,10)$ at time $t=5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the trajectory to verify that the solution is correct\n",
    "plt.xlim(-5, 105)\n",
    "plt.ylim(0, 40)\n",
    "plt.plot(X,Y,linewidth=4) # Plot the trajectory\n",
    "plt.plot([0,j],[0,h],'or',markersize=15) # Plot the start and end points as red circles\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we are going to approximate the solution $\\vec x(t)$ at discrete points in time.  In this block of code, you will need to fill in the entries for the matrix $A$ and the vector $\\vec b$.  See the problem description for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A=np.zeros((2*t.shape[0],2*t.shape[0]))\n",
    "b=np.ones((2*t.shape[0],1))  # Set up as b[2*i]=x_i, b[2*i+1]=y_i\n",
    "\n",
    "# Initial conditions\n",
    "# BEGIN STUDENT: Include the values in the vector b that correspond with the equations for the boundary conditions\n",
    "b[0,0]=\n",
    "b[1,0]=\n",
    "b[-2,0]=\n",
    "b[-1,0]=\n",
    "# END STUDENT\n",
    "\n",
    "# Set up the equations for the initial conditions in the A matrix\n",
    "# BEGIN STUDENT: Include the values in the matrix A that correspond with the equations for the boundary conditions\n",
    "A[0,0]=\n",
    "A[1,1]=\n",
    "A[-2,-2]=\n",
    "A[-1,-1]=\n",
    "# END STUDENT\n",
    "\n",
    "\n",
    "# Set up the acceleration approximation equations\n",
    "for i in range(1,t.shape[0]-1):\n",
    "    # BEGIN STUDENT: Include the values in the matrix A that correspond with the acceleration approximation equations\n",
    "    # These three entries correspond with the x component\n",
    "    A[2*i,2*i-2]   = \n",
    "    A[2*i,2*i+2]   = \n",
    "    A[2*i,2*i]     = \n",
    "    \n",
    "    # These three entries correspond with the y component\n",
    "    A[2*i+1,2*i-1] = \n",
    "    A[2*i+1,2*i+3] = \n",
    "    A[2*i+1,2*i+1] = \n",
    "    \n",
    "    # STUDENT: Include the values in the vector b for the acceleration approximation equations\n",
    "    b[2*i]   = \n",
    "    b[2*i+1] = \n",
    "    # END STUDENT\n",
    "    \n",
    "# Solve for the points\n",
    "Ainv = np.linalg.inv(A)\n",
    "p = np.asmatrix(Ainv) * np.asmatrix(b)\n",
    "\n",
    "# Get the answer\n",
    "Xapprox = np.zeros(X.shape)\n",
    "Yapprox = np.zeros(Y.shape)\n",
    "Xapprox.shape\n",
    "for i in range(0,t.shape[0]):\n",
    "    Xapprox[i] = p[2*i,0]\n",
    "    Yapprox[i] = p[2*i+1,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot our approximation as a thick blue line and plot the analytic solution as a thin magenta line on top.  If $A$ and $\\vec b$ are filled in correctly, then the two trajectories should overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the trajectory to verify that the solution is correct\n",
    "plt.xlim(-5, 105)\n",
    "plt.ylim(0, 40)\n",
    "plt.plot(Xapprox,Yapprox,linewidth=4) # Plot the trajectory\n",
    "plt.plot(X,Y,'-m',linewidth=2) # Plot the exact solution on top to compare\n",
    "plt.plot([0,j],[0,h],'or',markersize=15) # Plot the start and end points as red circles\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now estimate the initial velocity in the approximated solution.  The velocity can be estimated by $\\vec v(t_0) \\approx \\frac{\\vec x(t_1) - \\vec x(t_0)}{h}$.  Because this is just an approximation of the velocity, we don't expect this to be exactly equal to the analytic solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vxApprox = (Xapprox[1]-Xapprox[0])/step\n",
    "vyApprox = (Yapprox[1]-Yapprox[0])/step\n",
    "print('Actual vx:'+str(vx)+' vy:'+str(vy))\n",
    "print('Approximation vx:'+str(vxApprox)+' vy:'+str(vyApprox))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the new function of acceleration that depends on time.  Notice that the direction of gravity changes over time as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(t):\n",
    "    return g+2.0*np.exp(t/3.0)+8*np.sin(t*2.0)\n",
    "\n",
    "G=f(t)\n",
    "plt.plot(t,G,linewidth=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Update the vector b\n",
    "for i in range(1,t.shape[0]-1):\n",
    "    # BEGIN STUDENT: Fill in the new values in b to account for the new acceleration\n",
    "    b[2*i+1] = \n",
    "    # END STUDENT\n",
    "\n",
    "#Find the new points\n",
    "p = np.asmatrix(Ainv) * np.asmatrix(b)\n",
    "\n",
    "# Get the answer\n",
    "Xapprox = np.zeros(X.shape)\n",
    "Yapprox = np.zeros(Y.shape)\n",
    "Xapprox.shape\n",
    "for i in range(0,t.shape[0]):\n",
    "    Xapprox[i] = p[2*i,0]\n",
    "    Yapprox[i] = p[2*i+1,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we plot the new trajectory.  Note that without constant acceleration, the plot is no longer parabolic.  Like the previous plot, you can verify that your answer is correct if the trajectory starts and ends on the two red circles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the trajectory to verify that the solution is correct\n",
    "plt.xlim(-5, 105)\n",
    "plt.ylim(0, 40)\n",
    "plt.plot(Xapprox,Yapprox,linewidth=4) # Plot the trajectory\n",
    "plt.plot([0,j],[0,h],'or',markersize=15) # Plot the start and end points as red circles\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will give an approximation for the initial velocity of the projectile with this new trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vxApprox = (Xapprox[1]-Xapprox[0])/step\n",
    "vyApprox = (Yapprox[1]-Yapprox[0])/step\n",
    "print('Approximation vx:'+str(vxApprox)+' vy:'+str(vyApprox))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}