# MatrixExp.py: Calculate the exponential of a matrix.
#
import numpy as N
import numpy.linalg as LA
# Calculate exponential of a matrix, using its eigensystem
def matrixExponentialEigen(matrix):
eigenvalues, eigenvectors = LA.eig(matrix)
eigenvalues = N.exp(eigenvalues)
L = eigenvalues * N.identity(eigenvalues.shape[0])
return N.dot(N.dot(eigenvectors,L),LA.inv(eigenvectors))
# Calculate an exponential by Taylor expansion
def matrixExponentialTaylor(matrix,order):
sum = N.identity(matrix.shape[0])
product = sum
for i in range(1, order):
product = N.dot(product, matrix)/i
sum = sum + product
return sum
# The examples in the Homework:
# Uncomment the one you want.
#A = N.array( [ [ 2., 1.] , [1., 1.] ] )
#A = N.array([ [ 1., 0., 0.],[0.,2.,0.],[0.,0.,3.] ] )
A = N.array( [ [ 1., 2., 3.] , [ 3., 2., 1.] , [4., 5., 6.] ] )
# Compare the results
print "A:"
print A
print "e^A: Eigen method"
print matrixExponentialEigen(A)
print "e^A: Taylor"
print matrixExponentialTaylor(A,5)
Functions of a matrix: A very slight modification let's you
calculate any function of a matrix.
As a sanity check, we also include a range of test cases and then, finally, tests for the identity relation: e^A * e^-A = I.
# FunctionOfMatrix.py: How to calculate arbitrary function of a matrix
# Note: This is a direct generalization of MatrixExp.py that allows
# one to pass in any function. Otherwise, it is the same.
#
import numpy as N
import numpy.linalg as LA
# Calculate an arbitrary function of a matrix
def matrixFunction(function, matrix):
eigenvalues, eigenvectors = LA.eig(matrix)
eigenvalues = function(eigenvalues)
L = eigenvalues * N.identity(eigenvalues.shape[0])
return N.dot(N.dot(eigenvectors,L),LA.inv(eigenvectors))
# Calculate an exponential by Taylor expansion
def matrixExponential(matrix,order):
sum = N.identity(matrix.shape[0])
product = sum
for i in range(1, order):
product = N.dot(product, matrix)/i
sum = sum + product
return sum
# Create a symmetric matrix.
#A = N.reshape(N.arange(9.), (3, 3))**2/10.
#A = (A + N.transpose(A))/2.
# Some simple cases to test:
#A = N.identity(2)
#A = N.ones((2,2))
#A = N.array( [ [ 1., 1e-10] , [1e-10, 1.] ] )
#A = N.identity(3)
#A = N.ones((3,3))
# The examples in the Homework:
#A = N.array( [ [ 2., 1.] , [1., 1.] ] )
#A = N.array([ [ 1., 0., 0.],[0.,2.,0.],[0.,0.,3.] ] )
A = N.array( [ [ 1., 2., 3.] , [ 3., 2., 1.] , [4., 5., 6.] ] )
# Compare the results
print "A:"
print A
print "e^A: Eigen method"
print matrixFunction(N.exp, A)
print "e^A: Taylor"
print matrixExponential(A,5)
# Check e^A x e^-A = I
print "Inverse: Eigen method"
eA = matrixFunction(N.exp, A)
emA = matrixFunction(N.exp, -A)
print "-A ="
print -A
print "e^-A ="
print emA
print "e^A x e^-A = I?"
print N.dot(eA,emA)
# Taylor
print "Inverse: Taylor"
print "e^-A ="
eATay = matrixExponential(A,10)
emATay = matrixExponential(-A,10)
print emATay
print "e^A x e^-A = I?"
print N.dot(eATay,emATay)
# PlotMaxSep.py: Plot the Logistic map maximum separation time series
# computed by MaxSeparation.py.
# The data was stored in MaxSep.txt
import matrixio
from numpy import *
from pylab import *
MaxSep = matrixio.readMatrix('MaxSep.txt')
xlist = zeros(len(MaxSep))
ylist = zeros(len(MaxSep))
for i in xrange(len(MaxSep)):
xlist[i] = MaxSep[i][0]
ylist[i] = MaxSep[i][1]
plot(xlist,ylist)
xlabel('Iterations')
ylabel('Maximum Separation')
title('Logistic map: Spreading of initially close states')
show()
# LogisticICSpread.py: Plot the mean and standard deviation of an ensemble
# of states as it is iterated by the Logistic map at r = 4.0.
import numpy as np
from pylab import *
def LogisticMap(x,r,nIts):
for i in xrange(nIts):
x = r * x * (1 - x)
return x
# Store our calculated results
MaxIts = 50
Time = np.zeros(MaxIts)
Mean = np.zeros(MaxIts)
StdDev = np.zeros(MaxIts)
# Set up the ensemble of initial conditions (states), as requested
nPts = 1000
IC = 0.3
Width = 1e-3
Ensemble = np.random.normal(IC,Width,nPts)
for i in xrange(MaxIts):
for j in xrange(nPts):
Ensemble[j] = LogisticMap(Ensemble[j],4.0,1)
Time[i] = i
Mean[i] = np.mean(Ensemble)
StdDev[i] = np.std(Ensemble)
plot(Time,Mean,'g-',label='Ensemble mean',linewidth=2)
plot(Time,StdDev,'r-',label='Ensemble standard deviation',linewidth=2)
xlabel('Iterations')
ylabel('Mean and standard deviation')
axis([0,MaxIts,0,1])
title('Logistic map: Mean and standard deviation of point cloud')
legend()
show()
# LogisticDensity.py: Plot the distribution (density, really) of a
# large number of Logistic map iterations.
import numpy as np
from pylab import *
def LogisticMap(x,r,nIts):
for i in xrange(nIts):
x = r * x * (1 - x)
return x
nTransients = 100
nIterations = 10000
# Control parameter
r = 4.0
# Store the iterates
Its = np.zeros(nIterations)
x = 0.3
# Iterate away transients
x = LogisticMap(x,r,nTransients)
# Collect iterates
for i in xrange(nIterations):
x = LogisticMap(x,r,1)
Its[i] = x
nBins = 20
hg = np.histogram(Its,nBins,(0.0,1.0))
# Turn the bin counts into probabilities
Density = hg[0] / float(nIterations)
# Make the bin boundary array and the count array the same length
Density = np.insert(Density,0,Density[0])
plot(hg[1],Density,'g-',label='Density',linewidth=2)
xlabel('x')
ylabel('Prob(x)')
# Rescale vertical axis so that we can see some detail
MaxProb = np.max(Density)
axis([0.0,1.0,0.0,2.0*MaxProb])
title('Logistic map: Density over %d iterations' % (nIterations))
legend()
show()
# DensityEvolution.py: Plot the distribution (density, really) of an ensemble
# of states iterated by the Logistic map at r = 4.0.
import numpy as np
from pylab import *
def LogisticMap(x,r,nIts):
for i in xrange(nIts):
x = r * x * (1 - x)
return x
# The time to stop
nIterations = 10
# Control parameter
r = 4.0
# Set up the ensemble of initial conditions (states), as before
nPts = 10000
IC = 0.3
Width = 1e-3
Ensemble = np.random.normal(IC,Width,nPts)
for j in xrange(nPts):
Ensemble[j] = LogisticMap(Ensemble[j],r,nIterations)
nBins = 30
hg = np.histogram(Ensemble,nBins,(0.0,1.0))
# Turn the bin counts into probabilities
Density = hg[0] / float(nPts)
# Make the bin boundary array and the count arrays the same length
Density = np.insert(Density,0,Density[0])
plot(hg[1],Density,'g-',label='Density',linewidth=2)
xlabel('x')
ylabel('Prob(x)')
# Rescale vertical scale so that we can see some detail
MaxProb = np.max(Density)
axis([0.0,1.0,0.0,2.0*MaxProb])
title('Logistic map: Density (%d states) after %d iterations' % (nPts,nIterations))
legend()
show()