## module newtonPoly
''' p = evalPoly(a,xData,x).
Evaluates Newton's polynomial p at x. The coefficient
vector 'a' can be computed by the function 'coeffts'.
a = coeffts(xData,yData).
Computes the coefficients of Newton's polynomial.
'''
def evalPoly(a,xData,x):
n = len(xData) - 1 # Degree of polynomial
p = a[n]
for k in range(1,n+1):
p = a[n-k] + (x -xData[n-k])*p
return p
def coeffts(xData,yData):
m = len(xData) # Number of data points
a = yData.copy()
for k in range(1,m):
a[k:m] = (a[k:m] - a[k-1])/(xData[k:m] - xData[k-1])
return a
## module neville
''' p = neville(xData,yData,x).
Evaluates the polynomial interpolant p(x) that passes
trough the specified data points by Neville's method.
'''
def neville(xData,yData,x):
m = len(xData) # number of data points
y = yData.copy()
for k in range(1,m):
y[0:m-k] = ((x - xData[k:m])*y[0:m-k] + \
(xData[0:m-k] - x)*y[1:m-k+1])/ \
(xData[0:m-k] - xData[k:m])
return y[0]
## example3_4
# The data points in the table lie on the plot of f (x) = 4.8cosπx/20. Interpolate this
# data by Newton’s method at x = 0, 0.5, 1.0, . . . , 8.0, and compare the results with the
# “exact” values yi = f (xi ).
# x 0.15 2.30 3.15 4.85 6.25 7.95
# y 4.79867 4.49013 4.2243 3.47313 2.66674 1.51909
import numpy as np
import math
from newtonPoly import *
xData = np.array([0.15,2.3,3.15,4.85,6.25,7.95])
yData = np.array([4.79867,4.49013,4.2243,3.47313,2.66674,1.51909])
a = coeffts(xData,yData)
print(" x yInterp yExact")
print("-----------------------")
for x in np.arange(0.0,8.1,0.5):
y = evalPoly(a,xData,x)
yExact = 4.8*math.cos(math.pi*x/20.0)
print(’{:3.1f} {:9.5f} {:9.5f}’.format(x,y,yExact))
input("\nPress return to exit")
## example 3_6
# Interpolate the data shown at x-increments of 0.05 and plot the results. Use both the
# polynomial interpolation and the rational function interpolation.
# x 0.1 0.2 0.5 0.6 0.8 1.2 1.5
# y −1.5342 −1.0811 −0.4445 −0.3085 −0.0868 0.2281 0.3824
import numpy as np
from rational import *
from neville import *
import matplotlib.pyplot as plt
xData = np.array([0.1,0.2,0.5,0.6,0.8,1.2,1.5])
yData = np.array([-1.5342,-1.0811,-0.4445,-0.3085, \
-0.0868,0.2281,0.3824])
x = np.arange(0.1,1.55,0.05)
n = len(x)
y = np.zeros((n,2))
for i in range(n):
y[i,0] = rational(xData,yData,x[i])
y[i,1] = neville(xData,yData,x[i])
plt.plot(xData,yData,’o’,x,y[:,0],’-’,x,y[:,1],’--’)
plt.xlabel(’x’)
plt.legend((’Data’,’Rational’,’Neville’),loc = 0)
plt.show()
input("\nPress return to exit")
## module cubicSpline
''' k = curvatures(xData,yData).
Returns the curvatures of cubic spline at its knots.
y = evalSpline(xData,yData,k,x).
Evaluates cubic spline at x. The curvatures k can be
computed with the function 'curvatures'.
'''
import numpy as np
from LUdecomp3 import *
def curvatures(xData,yData):
n = len(xData) - 1
c = np.zeros(n)
d = np.ones(n+1)
e = np.zeros(n)
k = np.zeros(n+1)
c[0:n-1] = xData[0:n-1] - xData[1:n]
d[1:n] = 2.0*(xData[0:n-1] - xData[2:n+1])
e[1:n] = xData[1:n] - xData[2:n+1]
k[1:n] =6.0*(yData[0:n-1] - yData[1:n]) \
/(xData[0:n-1] - xData[1:n]) \
-6.0*(yData[1:n] - yData[2:n+1]) \
/(xData[1:n] - xData[2:n+1])
LUdecomp3(c,d,e)
LUsolve3(c,d,e,k)
return k
def evalSpline(xData,yData,k,x):
def findSegment(xData,x):
iLeft = 0
iRight = len(xData)- 1
while 1:
if (iRight-iLeft) <= 1: return iLeft
i =(iLeft + iRight)/2
if x <= xData[i]: iRight = i
else: iLeft = i
i = findSegment(xData,x)
h = xData[i] - xData[i+1]
y = ((x - xData[i+1])**3/h - (x - xData[i+1])*h)*k[i]/6.0 \
- ((x - xData[i])**3/h - (x - xData[i])*h)*k[i+1]/6.0 \
+ (yData[i]*(x - xData[i+1]) \
- yData[i+1]*(x - xData[i]))/h
return y
## example3_9
# program that interpolates, using the module cubicSpline, between given
# data points with a natural cubic spline. The program evaluate the interpolant
# for more than one value of x. Use data points specified f (x) = 4.8cos(πx/20)
# at x = 0, 0.5, 1.0, . . . , 8.0 and compute the interpolant at
# x = 1.5 and x = 4.5 (because of symmetry, these values should be equal)
import numpy as np
from cubicSpline import *
xData = np.array([1,2,3,4,5],float)
yData = np.array([0,1,0,1,0],float)
k = curvatures(xData,yData)
while True:
try: x = eval(input("\nx ==> "))
except SyntaxError: break
print("y =",evalSpline(xData,yData,k,x))
input("Done. Press return to exit")
## module polyFit
''' c = polyFit(xData,yData,m).
Returns coefficients of the polynomial
p(x) = c[0] + c[1]x + c[2]x^2 +...+ c[m]x^m
that fits the specified data in the least
squares sense.
sigma = stdDev(c,xData,yData).
Computes the std. deviation between p(x)
and the data.
'''
import numpy as np
import math
from gaussPivot import *
def polyFit(xData,yData,m):
a = np.zeros((m+1,m+1))
b = np.zeros(m+1)
s = np.zeros(2*m+1)
for i in range(len(xData)):
temp = yData[i]
for j in range(m+1):
b[j] = b[j] + temp
temp = temp*xData[i]
temp = 1.0
for j in range(2*m+1):
s[j] = s[j] + temp
temp = temp*xData[i]
for i in range(m+1):
for j in range(m+1):
a[i,j] = s[i+j]
return gaussPivot(a,b)
def stdDev(c,xData,yData):
def evalPoly(c,x):
m = len(c) - 1
p = c[m]
for j in range(m):
p = p*x + c[m-j-1]
return p
n = len(xData) - 1
m = len(c) - 1
sigma = 0.0
for i in range(n+1):
p = evalPoly(c,xData[i])
sigma = sigma + (yData[i] - p)**2
sigma = math.sqrt(sigma/(n - m))
return sigma
## module plotPoly
''' plotPoly(xData,yData,coeff,xlab='x',ylab='y')
Plots data points and the fitting
polynomial defined by its coefficient
array coeff = [a0, a1. ...]
xlab and ylab are optional axis labels
'''
import numpy as np
import matplotlib.pyplot as plt
def plotPoly(xData,yData,coeff,xlab='x',ylab='y'):
m = len(coeff)
x1 = min(xData)
x2 = max(xData)
dx = (x2 - x1)/20.0
x = np.arange(x1,x2 + dx/10.0,dx)
y = np.zeros((len(x)))*1.0
for i in range(m):
y = y + coeff[i]*x**i
plt.plot(xData,yData,'o',x,y,'-')
plt.xlabel(xlab); plt.ylabel(ylab)
plt.grid (True)
plt.show()
## example3_12
# program that fits a polynomial of arbitrary degree mto the data points shown
# in the following table. Use the program to determine m that best fits this data in the
# least-squares sense.
# x −0.04 0.93 1.95 2.90 3.83 5.00
# y −8.66 −6.44 −4.36 −3.27 −0.88 0.87
# x 5.98 7.05 8.21 9.08 10.09
# y 3.31 4.63 6.19 7.40 8.85
import numpy as np
from polyFit import *
xData = np.array([-0.04,0.93,1.95,2.90,3.83,5.0, \
5.98,7.05,8.21,9.08,10.09])
yData = np.array([-8.66,-6.44,-4.36,-3.27,-0.88,0.87, \
3.31,4.63,6.19,7.4,8.85])
while True:
try:
m = eval(input("\nDegree of polynomial ==> "))
coeff = polyFit(xData,yData,m)
print("Coefficients are:\n",coeff)
print("Std. deviation =",stdDev(coeff,xData,yData))
except SyntaxError: break
input("Finished. Press return to exit")
## example2_13
# function that inverts a matrix using LU decomposition with pivoting. Test the
# function by inverting
# A =[0.6 −0.4 1.0
# −0.3 0.2 0.5
# 0.6 −1.0 0.5]
import numpy as np
from LUpivot import *
def matInv(a):
n = len(a[0])
aInv = np.identity(n)
a,seq = LUdecomp(a)
for i in range(n):
aInv[:,i] = LUsolve(a,aInv[:,i],seq)
return aInv
a = np.array([[ 0.6, -0.4, 1.0],\
[-0.3, 0.2, 0.5],\
[ 0.6, -1.0, 0.5]])
aOrig = a.copy() # Save original [a]
aInv = matInv(a) # Invert [a] (original [a] is destroyed)
print("\naInv =\n",aInv)
print("\nCheck: a*aInv =\n", np.dot(aOrig,aInv))
input("\nPress return to exit")
## example2_14
# Invert thematrix
# A =[2 −1 0 0 0 0
# −1 2 −1 0 0 0
# 0 −1 2 −1 0 0
# 0 0 −1 2 −1 0
# 0 0 0 −1 2 −1
# 0 0 0 0 −1 5]
import numpy as np
from LUdecomp3 import *
n = 6
d = np.ones((n))*2.0
e = np.ones((n-1))*(-1.0)
c = e.copy()
d[n-1] = 5.0
aInv = np.identity(n)
c,d,e = LUdecomp3(c,d,e)
for i in range(n):
aInv[:,i] = LUsolve3(c,d,e,aInv[:,i])
print("\nThe inverse matrix is:\n",aInv)
input("\nPress return to exit")
## module gaussSeidel
''' x,numIter,omega = gaussSeidel(iterEqs,x,tol = 1.0e-9)
Gauss-Seidel method for solving [A]{x} = {b}.
The matrix [A] should be sparse. User must supply the
function iterEqs(x,omega) that returns the improved {x},
given the current {x} ('omega' is the relaxation factor).
'''
import numpy as np
import math
def gaussSeidel(iterEqs,x,tol = 1.0e-9):
omega = 1.0
k = 10
p = 1
for i in range(1,501):
xOld = x.copy()
x = iterEqs(x,omega)
dx = math.sqrt(np.dot(x-xOld,x-xOld))
if dx < tol: return x,i,omega
# Compute relaxation factor after k+p iterations
if i == k: dx1 = dx
if i == k + p:
dx2 = dx
omega = 2.0/(1.0 + math.sqrt(1.0 - (dx2/dx1)**(1.0/p)))
print('Gauss-Seidel failed to converge')
## module conjGrad
''' x, numIter = conjGrad(Av,x,b,tol=1.0e-9)
Conjugate gradient method for solving [A]{x} = {b}.
The matrix [A] should be sparse. User must supply
the function Av(v) that returns the vector [A]{v}
and the starting vector x.
'''
import numpy as np
import math
def conjGrad(Av,x,b,tol=1.0e-9):
n = len(b)
r = b - Av(x)
s = r.copy()
for i in range(n):
u = Av(s)
alpha = np.dot(s,r)/np.dot(s,u)
x = x + alpha*s
r = b - Av(x)
if(math.sqrt(np.dot(r,r))) < tol:
break
else:
beta = -np.dot(r,u)/np.dot(s,u)
s = r + beta*s
return x,i
## example2_17
# program to solve the following n simultaneous equations by the
# Gauss-Seidel method with relaxation (the program should work with any value of n):
# [2 −1 0 0 . . . 0 0 0 1
# −1 2 −1 0 . . . 0 0 0 0
# 0 −1 2 −1 . . . 0 0 0 0
# . . . . . . . .
# . . . . . . . .
# . . . . . . . .
# 0 0 0 0 . . . −1 2 −1 0
# 0 0 0 0 . . . 0 −1 2 −1
# 1 0 0 0 . . . 0 0 −1 2]
# [x1
# x2
# x3
# .
# .
# .
# xn−2
# xn−1
# xn] =
# [0
# 0
# 0
# .
# .
# .
# 0
# 0
# 1]
#Run the program with n = 20
import numpy as np
from gaussSeidel import *
def iterEqs(x,omega):
n = len(x)
x[0] = omega*(x[1] - x[n-1])/2.0 + (1.0 - omega)*x[0]
for i in range(1,n-1):
x[i] = omega*(x[i-1] + x[i+1])/2.0 + (1.0 - omega)*x[i]
x[n-1] = omega*(1.0 - x[0] + x[n-2])/2.0 \
+ (1.0 - omega)*x[n-1]
return x
n = eval(input("Number of equations ==> "))
x = np.zeros(n)
x,numIter,omega = gaussSeidel(iterEqs,x)
print("\nNumber of iterations =",numIter)
print("\nRelaxation factor =",omega)
print("\nThe solution is:\n",x)
input("\nPress return to exit")
## example2_18
# Solve Example 2.17 with the conjugate gradient method, also using n = 20
import numpy as np
from conjGrad import *
def Ax(v):
n = len(v)
Ax = np.zeros(n)
Ax[0] = 2.0*v[0] - v[1]+v[n-1]
Ax[1:n-1] = -v[0:n-2] + 2.0*v[1:n-1] -v[2:n]
Ax[n-1] = -v[n-2] + 2.0*v[n-1] + v[0]
return Ax
n = eval(input("Number of equations ==> "))
b = np.zeros(n)
b[n-1] = 1.0
x = np.zeros(n)
x,numIter = conjGrad(Ax,x,b)
print("\nThe solution is:\n",x)
print("\nNumber of iterations =",numIter)
input("\nPress return to exit")