newton and quasi-newton convergence experiments

Author: dreamchef

.DS_Store

@@ -0,0 +1,109 @@
+import numpy as np
+import numpy.linalg as la
+import matplotlib.pyplot as plt
+
+
+def driver():
+
+
+    f = lambda x: np.exp(x)
+
+    N = 3
+    ''' interval'''
+    a = 0
+    b = 1
+   
+   
+    ''' create equispaced interpolation nodes'''
+    xint = np.linspace(a,b,N+1)
+    
+    ''' create interpolation data'''
+    yint = f(xint)
+    
+    ''' create points for evaluating the Lagrange interpolating polynomial'''
+    Neval = 1000
+    xeval = np.linspace(a,b,Neval+1)
+    yeval_l= np.zeros(Neval+1)
+    yeval_dd = np.zeros(Neval+1)
+  
+    '''Initialize and populate the first columns of the 
+     divided difference matrix. We will pass the x vector'''
+    y = np.zeros( (N+1, N+1) )
+     
+    for j in range(N+1):
+       y[j][0]  = yint[j]
+
+    y = dividedDiffTable(xint, y, N+1)
+    ''' evaluate lagrange poly '''
+    for kk in range(Neval+1):
+       yeval_l[kk] = eval_lagrange(xeval[kk],xint,yint,N)
+       yeval_dd[kk] = evalDDpoly(xeval[kk],xint,y,N)
+          
+
+    
+
+
+    ''' create vector with exact values'''
+    fex = f(xeval)
+       
+
+    plt.figure()    
+    plt.plot(xeval,fex,'ro-')
+    plt.plot(xeval,yeval_l,'bs--') 
+    plt.plot(xeval,yeval_dd,'c.--')
+    plt.legend()
+
+    plt.figure() 
+    err_l = abs(yeval_l-fex)
+    err_dd = abs(yeval_dd-fex)
+    plt.semilogy(xeval,err_l,'ro--',label='lagrange')
+    plt.semilogy(xeval,err_dd,'bs--',label='Newton DD')
+    plt.legend()
+    plt.show()
+
+def eval_lagrange(xeval,xint,yint,N):
+
+    lj = np.ones(N+1)
+    
+    for count in range(N+1):
+       for jj in range(N+1):
+           if (jj != count):
+              lj[count] = lj[count]*(xeval - xint[jj])/(xint[count]-xint[jj])
+
+    yeval = 0.
+    
+    for jj in range(N+1):
+       yeval = yeval + yint[jj]*lj[jj]
+  
+    return(yeval)
+  
+    
+
+
+''' create divided difference matrix'''
+def dividedDiffTable(x, y, n):
+ 
+    for i in range(1, n):
+        for j in range(n - i):
+            y[j][i] = ((y[j][i - 1] - y[j + 1][i - 1]) /
+                                     (x[j] - x[i + j]));
+    return y;
+    
+def evalDDpoly(xval, xint,y,N):
+    ''' evaluate the polynomial terms'''
+    ptmp = np.zeros(N+1)
+    
+    ptmp[0] = 1.
+    for j in range(N):
+      ptmp[j+1] = ptmp[j]*(xval-xint[j])
+     
+    '''evaluate the divided difference polynomial'''
+    yeval = 0.
+    for j in range(N+1):
+       yeval = yeval + y[0][j]*ptmp[j]  
+
+    return yeval
+
+       
+
+driver()        

interp.py

@@ -0,0 +1,19 @@
+import numpy as np
+
+f = lambda x: x*x**2 - y**2
+g = lambda x: 3*x*y**2 - x**3 - 1
+
+x,y = 1,1
+
+matrix = [[1/6, 1/18],
+          [0, 1/6]]
+
+for i in range(1000):
+    vec = [x,y] - np.matmul(matrix, [f(x), g(y)])
+
+    x = vec[0]
+
+    y = vec[1]
+
+    print(x,y)
+

interp_system.py

@@ -0,0 +1,93 @@
+import numpy as np
+import math
+import time
+from numpy.linalg import inv
+from numpy.linalg import norm
+import matplotlib.pyplot as plt
+import sympy
+
+f = lambda x,y: x**2 + y**2 - 4
+g = lambda x,y: np.exp(x) + y - 1
+
+x,y = 1,1
+
+J = lambda x,y: [[2*x, 2*y],
+                    [np.exp(x), 1]]
+
+# J_broyden = lambda x,y:
+
+
+for init in [[1,1], [1,-1], [0,0]]:
+
+    print('\nINITIAL CONDIITONS', init)
+
+    x = init[0]
+    y = init[1]
+
+    # Newton's method
+
+    print('\nVANILLA NEWTON\n', x, y)
+
+
+    try:
+
+        for i in range(20):
+            vec = [x,y] - np.matmul(inv(J(x,y)), [f(x,y), g(x,y)])
+
+            x = vec[0]
+
+            y = vec[1]
+
+            print(x,y)
+
+    except:
+
+        print('Error')
+
+    
+
+
+    # Lazy Newton
+        
+    x = init[0]
+    y = init[1]
+        
+
+
+    print('\nLAZY\n', x, y)
+
+    # calc J once
+
+    J_ = J(x,y)
+
+
+    try:
+
+        for i in range(50):
+            vec = [x,y] - np.matmul(inv(J_), [f(x,y), g(x,y)])
+
+            x = vec[0]
+
+            y = vec[1]
+
+            print(x,y)
+
+    except:
+        
+        print('Error')
+        
+
+
+    x,y = init
+
+    print('\nBROYDEN\n', x, y)
+
+    for i in range(50):
+
+        vec = [x,y] - np.matmul(inv(J_broyden(x,y)), [f(x,y), g(x,y)])
+
+        x = vec[0]
+
+        y = vec[1]
+
+        print(x,y)