{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample Project Report: Find the equilbrium position of a charge\n",
    "### Author: Ilya Nemenman\n",
    "### Date: Feb 3, 2019\n",
    "### Other group members: None\n",
    "In this project, we solve the problem of finding the stationary position of a charge placed at an arbitrary point on a rod, with three other fixed charges of the same sign somewhere else on the rod. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem analysis\n",
    "### Physical considerations\n",
    "For the equilibrium position of the charge to exist, the charge must be between two other charges of the *same sign*, so that it gets repelled from both of them, and settles somewhere in the middle. The sign of the third fixed charge does not matter, as it will provide just a correction for the final position of the movng charge. There will be no equilibirum position if the moving charge is put to the left of the leftmost or the right of the rightmost fixed charge, and we check this in the code. Since the initial position of the moving charge can change depending on the input initial conditions, we will simplify the problem somewhat and will make an **assumption** (and check for it in the code) that all charges are of the same sign.\n",
    "\n",
    "An additional physical **assumption** is that there is friction, so that the charge actually settles down in its equilibroum position, rather than continues to oscillate forever.\n",
    "\n",
    "### Coding considerations\n",
    "The model being build is a *static* model for finding the equilibrium position of the charge. It's a deterministic model, since the charge equilibrium position is completely determined by the position and the value of the fixed charges. The variable of interest -- the equilibrium position is a continuous variable. To initialize the problem we will need to specify three positions of the fixed charges and the values of these charges. Additionally, we will need to specify the tolerance for finding the equilibrium soluiton. We will additionally count the number of iterations being performed, and if the change in the charge position after a certain number of iterations is still of order 1 or more, it means that the algorithm is not converging."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model development\n",
    "As discussed in class, we write down the second Newton's law for the moving charge, accounting for the three Coulomb forces from the stationary charges. We then set the force to zero -- the necessary condition for the equilibrium. Finally, after algebraic manipulations, we arrive at the following equiation for x, the equilibrium position of the charge:\n",
    "$$0=F(x)= q_1(x-p_2)^2(x-p_3)^2 - q_2(x-p_1)^2(x-p_3)^2 - q_2(x-p_2)^2(x-p_1)^2,$$\n",
    "where $q_i$ and $p_i$, $i=1,\\dots,3$ are the values and the position of the three fixed charges. This equation will be solved using the Newton-Raphson algorithm for root finding. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np #importing the numpy module\n",
    "\n",
    "p1 = 0 # position of three fixed charges\n",
    "p2 = 10 \n",
    "p3 = 12\n",
    "q1 = 1 # values of three charges\n",
    "q2 = 1\n",
    "q3 = 0\n",
    "x0 = 6 # initial condition\n",
    "Tolerance = 1e-5 # tolerance for the change in the charge position before termination\n",
    "MaxIter = 5 # number of iterations, after which we check for convergence "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Various exception checks for initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The moving charge position initial guess is between two other charges. A solution should exist.\n",
      "The fixed charges are all of the same sign. The solution should exist.\n"
     ]
    }
   ],
   "source": [
    "ErrVal = 0 # we will use ErrVal to keep track of whether there were any errors encountered.\n",
    "if np.min((p1,p2,p3))>=x0:\n",
    "    print('ERROR: The initial guess for the charge position is to the left of the leftmost fixed charge. The charge will run away to infinity.')\n",
    "    ErrVal = 1\n",
    "elif np.max((p1,p2,p3))<=x0:\n",
    "    print('ERROR: The initial guess for the charge position is to the right of the rightmost fixed charge. The charge will run away to infinity.')\n",
    "    ErrVal = 1\n",
    "else:\n",
    "    print('The moving charge position initial guess is between two other charges. A solution should exist.')\n",
    "if np.any(np.array((q1,q2))*q3<0): # this line checks if q3 is the same sign as both q1 and q2\n",
    "    print('ERROR: The fixed charges are not all the same sign. The solution may not exist.')\n",
    "    ErrVal = 1\n",
    "else:\n",
    "    print('The fixed charges are all of the same sign. The solution should exist.')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current change in the argument is -1.5.\n",
      "Current change in the argument is 0.4411764705882353.\n",
      "Current change in the argument is 0.05785920925747335.\n",
      "Current change in the argument is 0.0009640545745187745.\n",
      "Current change in the argument is 2.6557975214278305e-07.\n"
     ]
    }
   ],
   "source": [
    "x = x0 #initial value for starting the loop\n",
    "dx = 1.0 # fake \"last change in x\"\n",
    "n = 0 # number of iterations performed so far\n",
    "if(ErrVal==0): # only try to find a solution of no errors during initialization\n",
    "    while(np.abs(dx) > Tolerance): # note that we check for convergence based on the change \n",
    "        # of the charge position, not on the change of the function F that we are minimizing\n",
    "        n = n+1\n",
    "        F = q1*(x-p2)**2*(x-p3)**2 - q2*(x-p1)**2*(x-p3)**2 - q3*(x-p2)**2*(x-p1)**2\n",
    "        G = 2*q1*(x-p2)*(x-p3)**2 + 2*q1*(x-p2)**2*(x-p3) - 2*q2*(x-p1)*(x-p3)**2  -\\\n",
    "            2*q2*(x-p1)**2*(x-p3) - 2*q3*(x-p2)*(x-p1)**2 - 2*q3*(x-p2)**2*(x-p1)\n",
    "        dx = -F/G\n",
    "        x = x+dx\n",
    "        print(\"Current change in the argument is \", dx, \".\", sep=\"\")\n",
    "        if ((n>MaxIter)&(np.abs(dx)>0.01)):\n",
    "            print('ERROR: Maximum number of iterations reached, and no convergence detected. Exiting.')\n",
    "            ErrVal = 1\n",
    "            break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Termination and output of the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The charge will settle at 4.99999999999998.\n"
     ]
    }
   ],
   "source": [
    "if ErrVal:\n",
    "    print(\"Stationary position of the charge couldn't be found because of a previous error.\")\n",
    "else:\n",
    "    print(\"The charge will settle at \", x, \".\", sep=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Verification\n",
    "We have run the code above with different initial conditions, which has resulted in various exception checks.\n",
    "Additionally, we have run the code with the third charge being zero. This should result in the solution very similar to that of only two fixed charges. For example, $q_1=1$, $q_2=1$, and $q_3=0$ with $p_1=0$, $p_2$=10 results in $x=5$ by solving the quadratic equation and, with $x_0=6$ in $4.99999999999998$, as can be seen above, which is within the numerical accuracy.\n",
    "\n",
    "*Note that in the future, when we learn how to implement functions, you will be expected to provide verification by having a cell that calls various functions you defined to solve your model with different initial conditions, and plotting or otherwise visualizing the results. You will then be able to show results of verification as well as results of the main model solution at the same time.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "Solution for the equilibrium position of the charge was implemented. The program, based on the Newton-Raphson algorithm, finds the solution. Exceptions for position of the charges and the charge values were encoded. The program output has been verified against solution of a simpler problem with only two fixed charges, and the results matched.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}