small mods to Newton…

* added fdcoeffv.py to the repo as an extra python routine

Author: mspieg

09_ODE_ivp_part2.ipynb

@@ -2948,7 +2948,7 @@
     "\n",
     "$$\n",
     "\\begin{align}\n",
-    " \\mathbf{F}(\\mathbf{x}+\\boldsymbol{\\delta})&= \\mathbf{b} - A\\mathbf{x_0} - A\\boldsymbol{\\delta} = \\mathbf{0}\\\\\n",
+    " \\mathbf{F}(\\mathbf{x}_0+\\boldsymbol{\\delta})&= \\mathbf{b} - A\\mathbf{x_0} - A\\boldsymbol{\\delta} = \\mathbf{0}\\\\\n",
     "      &=\\mathbf{F}(\\mathbf{x}_0) - A\\boldsymbol{\\delta} = \\mathbf{0}\\\\\n",
     "\\end{align}\n",
     "$$"
@@ -3222,15 +3222,16 @@
    "outputs": [],
    "source": [
     "# Some quick and dirty Code\n",
-    "def newton(F,J,x0,tol=1.e-6,MAX_ITS=100):\n",
+    "def newton(F,J,x0,tol=1.e-6,MAX_ITS=100,verbose=True):\n",
     "    \"\"\" Solve F(x) = 0 using Newton's method until ||F(x)|| < tol or reaches Max iterations\"\"\"\n",
     "    \n",
     "    x = x0\n",
     "    for k in range(MAX_ITS+1):\n",
     "        delta = numpy.linalg.solve(J(x),-F(x))\n",
     "        x += delta\n",
     "        err = numpy.linalg.norm(F(x))\n",
-    "        print('||F|| = {}'.format(err))\n",
+    "        if verbose:\n",
+    "            print('||F|| = {}'.format(err))\n",
     "        if err < tol: \n",
     "            return x, k\n",
     "                \n",
@@ -3310,13 +3311,15 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {
     "hide_input": false
    },
-   "outputs": [],
-   "source": []
+   "source": [
+    "### An Exercise:  \n",
+    "\n",
+    "Write a BDF-2 scheme that can solve a general non-linear problem and test it using the vander-pol oscillator."
+   ]
   }
  ],
  "metadata": {

fdcoeffV.py

@@ -0,0 +1,64 @@
+import numpy as numpy
+from  scipy.special import factorial
+
+def fdcoeffV(k,xbar,x):
+    """
+    fdcoeffV routine modified from Leveque (2007) matlab function
+    
+    Params:
+    -------
+    
+    k: int
+        order of derivative
+    xbar: float
+        point at which derivative is to be evaluated
+    x: ndarray
+        numpy array of coordinates to use in calculating the weights
+    
+    Returns:
+    --------
+    c: ndarray
+        array of floats of coefficients.  
+
+    Compute coefficients for finite difference approximation for the
+    derivative of order k at xbar based on grid values at points in x.
+
+    WARNING: This approach is numerically unstable for large values of n since
+    the Vandermonde matrix is poorly conditioned.  Use fdcoeffF.m instead,
+    which is based on Fornberg's method.
+
+     This function returns a row vector c of dimension 1 by n, where n=length(x),
+     containing coefficients to approximate u^{(k)}(xbar), 
+     the k'th derivative of u evaluated at xbar,  based on n values
+     of u at x(1), x(2), ... x(n).  
+
+     If U is an array containing u(x) at these n points, then 
+     c.dot(U) will give the approximation to u^{(k)}(xbar).
+
+     Note for k=0 this can be used to evaluate the interpolating polynomial 
+     itself.
+
+    Requires len(x) > k.  
+    Usually the elements x(i) are monotonically increasing
+    and x(1) <= xbar <= x(n), but neither condition is required.
+    The x values need not be equally spaced but must be distinct.  
+    
+    Modified rom  http://www.amath.washington.edu/~rjl/fdmbook/  (2007)
+    """
+    
+
+    n = x.shape[0]
+    assert  k < n, " The order of the derivative must be less than the stencil width"
+
+    # Generate the Vandermonde matrix from the Taylor series
+    A = numpy.ones((n,n))
+    xrow = (x - xbar)  # displacements x-xbar 
+    for i in range(1,n):
+        A[i,:] = (xrow**(i))/factorial(i);
+        
+    b = numpy.zeros(n)    # b is right hand side,
+    b[k] = 1              # so k'th derivative term remains
+
+    c = numpy.linalg.solve(A,b)          # solve n by n system for coefficients
+    
+    return c