Miscellaneous cleanups.

2_num_lin_alg/chol_time.py

@@ -12,24 +12,24 @@
     a=np.random.random((m,m))
     b=np.dot(a,a.T)
 
-    # Store the initial wall clock time
-    t=time()
+    # Store the initial process time
+    t=process_time()
 
-    # Measure the process time
-    e=process_time()
+    # Measure the wall clock time
+    e=time()
     n=0
     f=e
 
     # Do as many Cholesky factorizations as possible within a 0.2 second
-    # interval of process time
+    # interval of wall clock time
     while f-e<0.2:
         n+=1
         l=scipy.linalg.cholesky(b)
-        f=process_time()
+        f=time()
 
-    # Print the average process time and wall clock time to do a Cholesky
+    # Print the average wall clock time and process time to do a Cholesky
     # factorization
-    print("%d %d %.5g %.5g"%(m,n,(f-e)/n,(time()-t)/n))
+    print("%d %d %.5g %.5g"%(m,n,(f-e)/n,(process_time()-t)/n))
 
     # Increase m by a multiplicative factor
     m+=m//10

2_num_lin_alg/lu_time.py

@@ -10,7 +10,7 @@
     # Construct a random matrix to apply the LU factorization to
     a=np.random.random((m,m))
 
-    # Store the initial wall clock time
+    # Store the initial process time
     t=process_time()
 
     # Measure the wall clock time

2_num_lin_alg/timing_test.py

@@ -1,21 +1,21 @@
 #!/usr/bin/python
-from time import process_time,time
+from time import time,process_time
 import numpy as np
 import scipy.linalg
 
-# Get the current time() (a measure of "wall-clock time", the physical passage of
-# time as you would observe by looking at a wall clock)
+# Get the current time(), a measure of "wall-clock time", the physical passage of
+# time as you would observe by looking at a wall clock
 t0=time()
 
-# Get the current process_time() (a measure of the time that this program has actually
-# spent on the processor)
+# Get the current process_time(), a measure of the time that this program has actually
+# spent on the processor
 p0=process_time()
 
 # Carry out a large, arbitrary calculation
 k=0
 for i in range(8000):
-	for j in range(8000):
-		k+=i+j
+    for j in range(8000):
+        k+=i+j
 
 # Print the elapsed time, using both the time() routine and the process_time() routine
 ttotal=time()-t0

3_num_calculus/diff.py

@@ -1,7 +1,6 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 import matplotlib.pyplot as plt
 import numpy as np
-import scipy.linalg
 from math import *
 
 # Function to consider
@@ -13,7 +12,7 @@ def df(z):
     return 10*sin(10*z)/(2+cos(10*z))**2
 
 # Differentiation matrix
-n=100
+n=200
 h=1/float((n-1))
 D=np.diag(-np.ones(n)/h)+np.diag(np.ones(n-1)/h,1)
 
@@ -28,15 +27,22 @@ def df(z):
 # Plot function
 x=np.linspace(0,1,n)
 y=np.array([1/(2+cos(10*xx)) for xx in x])
+plt.xlabel('x')
+plt.ylabel('f(x)')
 plt.plot(x,y)
 plt.show()
 
 # Calculate derivative and plot
 dy=np.dot(D,y)
-plt.plot(x,y,x,dy)
+plt.xlabel('x')
+plt.plot(x,y,label='f(x)')
+plt.plot(x,dy,label="f'(x)")
+plt.legend()
 plt.show()
 
 # Plot error
 err=np.array([dy[i] - df(x[i]) for i in range(n)])
+plt.xlabel('x')
+plt.ylabel('Derivative error')
 plt.plot(x,err)
 plt.show()

3_num_calculus/l-v.py

@@ -1,38 +1,38 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import odeint
-from math import *
 
 # Constants in model
-alpha=1.5
-beta=0.5
-gamma=0.4
-delta=0.4
+alpha1=1.2
+beta1=0.4
+alpha2=0.2
+beta2=0.1
 
 # Function that evaluates the RHS of the ODE. It has two components,
 # representing the changes in prey and predator populations.
-def deriv(x,t):
-    return np.array([alpha*x[0]-beta*x[0]*x[1],-gamma*x[1]+delta*x[0]*x[1]])
+def f(y,t):
+    return np.array([alpha1*y[0]-beta1*y[0]*y[1], \
+                     -alpha2*y[1]+beta2*y[0]*y[1]])
 
 # Specify the range of time values where the ODE should be solved at
 time=np.linspace(0,70,500)
 
 # Initial conditions, set to the initial populations of prey and predators
-xinit=np.array([10,5])
+yinit=np.array([10,5])
 
 # Solve ODE using the "odeint" library in SciPy
-x=odeint(deriv,xinit,time)
+y=odeint(f,yinit,time)
 
 # Print the solutions
-for i in range(0,500):
-    print((time[i],x[i,0],x[i,1]))
+#for i in range(0,500):
+#    print(time[i],x[i,0],x[i,1])
 
 # Plot the solutions
 plt.figure()
-p0,=plt.plot(time,x[:,0])
-p1,=plt.plot(time,x[:,1])
-plt.legend([p0,p1],["B(t)","C(t)"])
+p0,=plt.plot(time,y[:,0])
+p1,=plt.plot(time,y[:,1])
+plt.legend([p0,p1],["Prey","Predators"])
 plt.xlabel('t')
 plt.ylabel('Population')
 plt.show()

3_num_calculus/order2.py

@@ -1,11 +1,10 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 # Import math functions
 from math import *
 
 # Initial variables and constants
 t=0
 h=0.05
-lam=2
 
 # Function to integrate
 def f(t,y):

3_num_calculus/quadrat.py

@@ -1,28 +1,28 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from math import *
 
 # Function to integrate
 def f(x):
     return sin(x)
 
-# Trapezoidal rule
+# Composite trapezoid rule
 def trap(f,a,b,n):
 
-    # Step size
+    # Subinterval width
     h=(b-a)/float(n)
 
-    # Trapezoidal formula
+    # Trapezoid formula
     fi=0.5*(f(a)+f(b))
     for i in range(1,n):
         fi+=f(a+i*h)
 
     # Return scaled answer
     return fi*h
 
-# Simpson's rule
+# Composite Simpson rule
 def simp(f,a,b,n):
 
-    # Step size
+    # Subinterval width
     h=(b-a)/float(n)
 
     # Simpson's formula
@@ -37,14 +37,14 @@ def simp(f,a,b,n):
 # Exact answer
 ex=cos(1)-cos(5)
 
-# Integrate with grids of different sizes. The errors for the trapezoidal rule
-# decay like O(h^2), matching the error bound derivation presented in the
-# lecture. The errors for Simpson's rule decay like O(h^4), which is actually
-# better than the expected O(h^3) based on the derivation in the lecture. It
-# turns out the the h^3 error terms cancel, leading to better performance,
-# which can be proved with a more careful analysis. Note however that there is
-# nothing wrong with the O(h^3) bound derived in the lecture - all it gives is
-# an upper bound, which is still true.
+# Integrate with different numbers of subintervals. The errors for the
+# trapezoidal rule decay like O(h^2), matching the error bound derivation
+# presented in the lecture. The errors for Simpson's rule decay like O(h^4),
+# which is actually better than the expected O(h^3) based on the derivation in
+# the lecture. It turns out the the h^3 error terms cancel, leading to better
+# performance, which can be proved with a more careful analysis. Note however
+# that there is nothing wrong with the O(h^3) bound derived in the lecture -
+# all it gives is an upper bound, which is still true.
 j=1
 while j<=65536:
     print(j,4./float(j),trap(f,1,5,j)-ex,simp(f,1,5,j)-ex)

3_num_calculus/stiff.py

@@ -15,7 +15,7 @@
 # Starting time and timestep (currently chosen within the stability region of
 # the explicit method)
 t=0
-dt=0.001
+dt=0.1
 
 while t<2:
 

3_num_calculus/transp.py

@@ -11,7 +11,7 @@
 z=np.zeros((m,snaps+1))
 
 # PDE-related constants; try switching c to -0.1 to see the unstable scheme
-c=0.1
+c=-0.1
 dx=1.0/m
 dt=0.01
 nu=c*dt/dx

4_optimization/bisection.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from math import *
 
 # Function to consider

4_optimization/fsolve.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 import numpy as np
 from scipy.optimize import *
 

4_optimization/h_bfgs.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 import numpy as np
 import scipy.optimize as op
 

4_optimization/h_newton.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 import numpy as np
 
 # Himmelblau function - a common example used for testing optimization algorithms

4_optimization/iter.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from math import *
 
 def f(x):
@@ -9,6 +9,6 @@ def f(x):
 
 for i in range(20):
     xa=log(xa+2)
-    xb=exp(xb)-2
+    #xb=exp(xb)-2
 
-    print(xa,f(xa),xb,f(xb))
+    print(xa,f(xa))#,xb,f(xb))

4_optimization/iter_2d.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from math import sqrt
 
 def f(x1,x2):

4_optimization/linprog.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from cvxopt import matrix,solvers
 
 # Constraints

4_optimization/linprog_alt.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from scipy.optimize import linprog
 import numpy as np
 

4_optimization/nr_secant.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 from math import *
 
 def f(x):