Initial commit

Initial commit

Homework 2/HW2a_Problem1_cjf2123.py

@@ -0,0 +1,27 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Oct  2 22:12:56 2012
+
+@author: cjf2123
+"""
+#Chris Fong
+#APMA 4301 
+#Homework 2a Problem 1
+import scipy
+import sympi
+
+from fdcoeffV import fdcoeffV
+
+
+k = [1,2];
+h = 1;
+xbar = [0*h,0.5*h,1*h]
+x = [0*h,1*h,2*h]
+c = scipy.zeros((len(k)*len(xbar),len(x)))
+for i in range(len(k)):
+    for j in range(len(xbar)):
+        c[(i)*len(xbar)+(j),:] = fdcoeffV(k[i],xbar[j],x)
+        if i == 0:
+            print 'First derivative, xbar = %2.1f: ' %xbar[j] ,str(c[(i)*len(xbar)+(j),:]),'\r'
+        else:
+            print 'Second derivative, xbar = %2.1f: '

Homework 2/README.md

@@ -0,0 +1 @@
+# This is my README

Homework 2/cjf2123_APMA4301_HW2a_P2.py

@@ -0,0 +1,149 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Oct  3 17:08:19 2012
+
+@author: cjf2123
+"""
+#Chris Fong
+#APMA 4301 
+#Homework 2a Problem 1
+
+import scipy as sp
+import numpy as np
+import pylab
+import diffMatrix as dm
+
+
+
+#N = array([8, 16, 32, 64, 128, 256, 512, 1024])
+#M = array([8, 16, 32, 64, 128, 256, 512, 1024, 2048])
+N = array([128])
+M = array([128.0])
+stencil_size = array([3,5])
+h = 1./M
+eh13 = sp.zeros((len(N),len(h)))
+eh23 = sp.zeros((len(N),len(h)))
+eh15 = sp.zeros((len(N),len(h)))
+eh25 = sp.zeros((len(N),len(h)))
+ehQ1_13 = sp.zeros((len(N),len(h)))
+ehQ1_23 = sp.zeros((len(N),len(h)))
+for i in range(len(N)):
+
+    x = np.linspace(0,1,N[i]+1)
+    f = x**2 + np.sin(4*pi*x)
+    fp = 2*x + 4*(np.pi)*(np.cos(4*pi*x))
+    fpp = 2 - 16*(np.pi**2)*(np.sin(4*pi*x))
+    D1f3 = dm.setDn(1,x,3)*f
+    D2f3 = dm.setDn(2,x,3)*f
+    D1f5 = dm.setDn(1,x,5)*f
+    D2f5 = dm.setDn(2,x,5)*f
+
+    #compute error between problem 1 and 2  
+    fppp = -64*(np.pi**3)*(np.cos(4*pi*x))
+    for j in range(len(h)):
+    #compute absolute mesh error
+        eh13[i,j] = np.sqrt(h[j])*(np.linalg.norm(D1f3-fp))
+        eh23[i,j] = np.sqrt(h[j])*(np.linalg.norm(D2f3-fpp))
+        eh15[i,j] = np.sqrt(h[j])*(np.linalg.norm(D1f5-fp))
+        eh25[i,j] = np.sqrt(h[j])*(np.linalg.norm(D2f5-fpp))
+
+        #compute error between problem 1 and 2  
+        error_Q1_1D = (h[j]**3)/6*fppp      #First error term
+        error_Q1_2D = (h[j]**2)/12*fppp
+        ehQ1_13[i,j] = np.sqrt(h[j])*(np.linalg.norm(error_Q1_1D))
+        ehQ1_23[i,j] = np.sqrt(h[j])*(np.linalg.norm(error_Q1_2D))
+
+eh13 = np.transpose(eh13)         #N increases as down the rows
+eh23 = np.transpose(eh23)
+eh15 = np.transpose(eh15)         
+eh25 = np.transpose(eh25)  
+ehQ1_13 = np.transpose(ehQ1_13)         
+ehQ1_23 = np.transpose(ehQ1_23)
+
+#Problem 2 (a) (i)
+fig = pylab.figure()
+eh1x = fig.add_subplot(2,1,1)
+pylab.title("Problem 2_a_i: 1st derivative- Log-Log Plot of h vs eh")
+eh2x = fig.add_subplot(2,1,2)
+pylab.title("Problem 2_a_i: 2nd derivative: Log-Log Plot of h vs eh")
+line1 = eh1x.plot(h,eh13)
+line2 = eh2x.plot(h,eh23)
+eh1x.set_yscale('log')
+eh1x.set_xscale('log')
+eh2x.set_yscale('log')
+eh2x.set_xscale('log')
+eh1x.legend(["N=8","N=16","N=32","N=64","N=128","N=256","N=512","N=1024"],loc="best")
+eh2x.legend(["N=8","N=16","N=32","N=64","N=128","N=256","N=512","N=1024"],loc="best")
+eh1x.set_xlabel("Step Size (h)")
+eh2x.set_xlabel("Step Size (h)")
+eh1x.set_ylabel("Error (eh)")
+eh2x.set_ylabel("Error (eh)")
+pylab.show()
+
+#Problem 2 (a) (ii)
+hlog = np.log10(h)
+eh13log = np.log10(eh13)
+eh23log = np.log10(eh23)
+(p13,b13) = np.polyfit(hlog,eh13log,1)
+(p23,b23) = np.polyfit(hlog,eh23log,1)
+
+print 'Problem 2_a_ii: Order of convergence of the Error, 1st derivative, : ', str(p13), ' for h = ', str(h),'respectively. \r'
+print 'Problem 2_a_ii: Order of convergence of the Error, 2nd derivative, : ', str(p23), ' for h = ', str(h),'respectively. \r'
+
+
+#Problem 2 (a) (iii)
+#The first error term for the centered difference we found from problem 1 was
+#h^3/6*u'''(x) for the first derivative and 
+#h^2/12*u'''(x) for the second derivative
+#The 2-norm of this vector will be compared with our calculated setD*f - f',f''
+fig = pylab.figure()
+eh1x = fig.add_subplot(2,1,1)
+pylab.title("Problem 2_a_iii: 1st derivative- Log-Log Plot of h vs eh")
+eh2x = fig.add_subplot(2,1,2)
+pylab.title("Problem 2_a_iii: 2nd derivative: Log-Log Plot of h vs eh")
+line1 = eh1x.plot(h,ehQ1_13)
+line2 = eh2x.plot(h,ehQ1_23)
+eh1x.set_yscale('log')
+eh1x.set_xscale('log')
+eh2x.set_yscale('log')
+eh2x.set_xscale('log')
+eh1x.legend(["N=8","N=16","N=32","N=64","N=128","N=256","N=512","N=1024"],loc="best")
+eh2x.legend(["N=8","N=16","N=32","N=64","N=128","N=256","N=512","N=1024"],loc="best")
+eh1x.set_xlabel("Step Size (h)")
+eh2x.set_xlabel("Step Size (h)")
+eh1x.set_ylabel("Error (eh)")
+eh2x.set_ylabel("Error (eh)")
+
+
+#Problem 2 (b) (i)
+fig = pylab.figure()
+eh1x = fig.add_subplot(2,1,1)
+pylab.title("Problem 2_b_i: 1st derivative- Log-Log Plot of h vs eh")
+eh2x = fig.add_subplot(2,1,2)
+pylab.title("Problem 2_b_i: 2nd derivative: Log-Log Plot of h vs eh")
+line1 = eh1x.plot(h,eh15)
+line2 = eh2x.plot(h,eh25)
+eh1x.set_yscale('log')
+eh1x.set_xscale('log')
+eh2x.set_yscale('log')
+eh2x.set_xscale('log')
+
+eh1x.legend(["N=8","N=16","N=32","N=64","N=128","N=256","N=512","N=1024"],loc="best")
+eh2x.legend(["N=8","N=16","N=32","N=64","N=128","N=256","N=512","N=1024"],loc="best")
+eh1x.set_xlabel("Step Size (h)")
+eh2x.set_xlabel("Step Size (h)")
+eh1x.set_ylabel("Error (eh)")
+eh2x.set_ylabel("Error (eh)")
+pylab.show()
+
+#Problem 2 (b) (ii)
+hlog = np.log10(h)
+eh15log = np.log10(eh15)
+eh25log = np.log10(eh25)
+(p15,b15) = np.polyfit(hlog,eh15log,1)
+(p25,b25) = np.polyfit(hlog,eh25log,1)
+
+print 'Problem 2_b_ii: Order of convergence of the Error, 1st derivative, : ', str(p15), ' for h = ', str(h),'respectively. \r'
+print 'Problem 2_b_ii: Order of convergence of the Error, 2nd derivative, : ', str(p25), ' for h = ', str(h),'respectively. \r'
+
+

Homework 2/diffMatrix.py

@@ -0,0 +1,131 @@
+# -*- coding: utf-8 -*-
+"""
+diffMatrix:  example program to set up sparse differentiation matrix
+Created on Tue Sep 18 01:27:13 2012
+
+@author: mspieg
+"""
+
+import scipy
+import scipy.sparse as sp
+import numpy as np
+import pylab
+
+from fdcoeffF import fdcoeffF
+
+def setD(k,x):
+    """ 
+    example function for setting k'th order sparse differentiation matrix over 
+    arbitrary mesh x with a 3 point stencil
+    
+    input:
+        k = degree of derivative <= 2
+        x = numpy array of coordinates >=3 in length
+    returns:
+        D sparse differention matric
+    """
+
+    assert(k < 3) # check to make sure k < 3
+    assert(len(x) > 2)
+
+    N = len(x)
+    # initialize a sparse NxN matrix in "lil" (linked list) format
+    #D = sp.lil_matrix((N,N))
+    Dx = scipy.zeros((N,N))
+    # assign the one-sided k'th derivative at x[0]
+    Dx[0,0:3] = fdcoeffF(k,x[0],x[0:3])
+    #D[0,0:3] = fdcoeffF(k,x[0],x[0:3])
+    # assign centered k'th ordered derivatives in the interior
+    for i in xrange(1,N-1):
+        #D[i,i-1:i+2] = fdcoeffF(k,x[i],x[i-1:i+2])
+        Dx[i,i-1:i+2] = fdcoeffF(k,x[i],x[i-1:i+2])
+    # assign one sided k'th derivative at end point x[-1]
+    #D[N-1,-3:] = fdcoeffF(k,x[N-1],x[-3:])
+    #Dx[N-1,-3:] = fdcoeffF(k,x[N-1],x[-3:])
+
+    # convert to csr (compressed row storage) and return
+    return D.tocsr()
+
+
+def setDn(k,x,ssize):
+    """ 
+    example function for setting k'th order sparse differentiation matrix over 
+    arbitrary mesh x with a 5 point stencil
+    
+    input:
+        k = degree of derivative <= 2
+        x = numpy array of coordinates >=5 in length
+    returns:
+        D sparse differention matric
+    """
+    N = len(x)
+    assert(k < 3) # check to make sure k < 3
+    assert(N > ssize+1)
+
+
+    mid = int(np.floor(ssize/2))
+    # initialize a sparse NxN matrix in "lil" (linked list) format
+    D = sp.lil_matrix((N,N))
+#    Dx = scipy.zeros((N,N))
+    # assign the one-sided k'th derivative at x[0]
+    for i in xrange(0,mid):
+        D[i,0:ssize] = fdcoeffF(k,x[i],x[i:ssize+i])
+#        Dx[i,0:ssize] = fdcoeffF(k,x[i],x[i:ssize+i])
+
+    # assign centered k'th ordered derivatives in the interior
+    for i in xrange(mid,N-mid):
+        D[i,i-mid:i+mid+1] = fdcoeffF(k,x[i],x[i-mid:i+mid+1])
+#        Dx[i,i-mid:i+mid+1] = fdcoeffF(k,x[i],x[i-mid:i+mid+1])
+
+    # assign one sided k'th derivative at end point x[-1]
+    for i in xrange(N-mid,N):
+        D[i,N-ssize:N] = fdcoeffF(k,x[i],x[i-ssize+1:i+1])
+#        Dx[i,N-ssize:N] = fdcoeffF(k,x[i],x[i-ssize+1:i+1])
+
+    # convert to csr (compressed row storage) and return
+    return D.tocsr()
+
+
+# quicky test program
+def plotDf(x,f,title=None):
+    """ Quick test routine to display derivative matrices and plot
+        derivatives of an arbitrary function f(x)
+        input: x: numpy array of mesh-points
+               f: function pointer to function to differentiate
+    """  
+    # calculate first and second derivative matrices
+    D1 = setDn(1,x,5) 
+#    D1 = setD(1,x) 
+    D2 = setD(2,x)
+#    
+#    print D1
+#    print D2
+
+    # show sparsity pattern of D1
+    pylab.figure()
+    pylab.spy(D1)
+    pylab.title("Sparsity pattern of D1")
+
+    # plot a function and it's derivatives
+    y = f(x)    
+    pylab.figure()
+    pylab.plot(x,y,x,D1*y,x,D2*y)
+    pylab.legend(['f','D1*f','D2*f'],loc="best")
+    if title:
+        pylab.title(title)
+    pylab.show()
+
+
+def main():
+    # set numpy grid array to be evenly spaced    
+    x = np.linspace(0,2,12)
+    # choose a simple quadratic function     
+    def f(x):
+        return x**2 + x
+
+    plotDf(x,f,"f=x^2 + x")
+
+
+
+if __name__ == "__main__":
+    main()

Homework 3/cjf2123_PS2b_Question1.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Oct 10 14:12:09 2012
+
+@author: cjf2123
+"""
+import scipy
+import scipy.sparse as sp
+import numpy as np
+#from scipy.sparse.linalg import spsolve
+#import pylab
+#from mpl_toolkits.mplot3d import Axes3D
+
+a = 0.0
+b = 1.0
+N = 100
+h = 1/(N-1.0)
+n = 2
+m = 2
+I = sp.eye(N,N)
+g = np.ones(N)
+T = sp.spdiags([g,-4*g,g],[-1,0,1],N,N)
+S = sp.spdiags([g,g],[-1,1],N,N)
+A = (sp.kron(I,T) + sp.kron(S,I)) / (h**2)
+A = A.tocsr()
+u = np.zeros((N,N))
+for i in range(N):
+    for j in range(N):
+        u[i,j] = np.sin(m*np.pi*i*h)*np.sin(n*np.pi*j*h)
+uvec = u.reshape((N*N,1))
+duvec = A*uvec
+du = duvec.reshape((N,N))   
+p=0
+
+