Chapter 14 - Boundary value and eigen value problems

Example No. 14_01 Pg No. 467

In [12]:
#Shooting Method

def heun(f,x0,y0,z0,h,xf):
    x = [x0] #
    global y
    y = [y0] #
    z = [z0] #
    n = (xf - x0)/h
    m1=[0,0];m2=[0,0];m=[0,0]
    for i in range(0,int(n)):
        m1[0] = z[i] 
        m1[1] = f(x[i],y[i])
        m2[0] = z[i] + h*m1[1]
        m2[1] = f(x[i]+h,y[i]+h*m1[0])
        m[0] = (m1[0] + m2[0])/2 
        m[1] = ( m1[1] + m2[1] )/2
        x.append(x[(i)] + h)
        y.append(y[(i)] + h*m[0])
        z.append(z[(i)] + h*m[1])
    
    B = y[int(n)]
    return B


def f(x,y):
    F = 6*x
    return F

x0 = 1 #
y0 = 2 #
h = 0.5 #
z0 = 2
M1 = z0 
xf = 2
B = 9
B1 = heun(f,x0,y0,z0,h,xf)
print 'B1 = ',B1

if B1 != B:
    print 'Since B1 is less than B , let z(1) = y(1) = 4*(M2)'
    z0 = 4
    M2 = z0
    B2 = heun(f,x0,y0,z0,h,xf)
    print 'B2 = ',B2
    if B2 != B:
        print 'Since B2 is larger than B ,let us have third estimate of z(1) = M3 '
        M3 = M2 - (B2 - B)*(M2 - M1)/(B2 - B1)
        z0 = M3 
        B3= heun(f,x0,y0,z0,h,xf)
        print 'B3 = ',B3
        print 'The solution is ',y

B1 =  7.75
Since B1 is less than B , let z(1) = y(1) = 4*(M2)
B2 =  9.75
Since B2 is larger than B ,let us have third estimate of z(1) = M3 
B3 =  9.0
The solution is  [2, 4.375, 9.0]

Example No. 14_02 Pg No. 470

In [1]:
from numpy import array,exp,zeros,vstack,hstack,linalg
#Finite Difference Method

def d2y(x):
    D2Y = exp(x**2)
    return D2Y
x_1 = 0#
y_0 = 0 #
y_1  = 0 #
h = 0.25
xf = 1
n = (xf-x_1)/h
A=zeros([int(n-1),int(n-1)])
B=zeros([int(n-1),1])
for i in range(0,int(n-1)):
    A[i,:] = [1,  -2,   1]
    B[i,0] = exp((x_1 + i*h)**2)*h**2

A[0,0] = 0 # #since we know y0 and y4
A[2,2] = 0 #
A[0,:] = hstack([ A[0,1:3], [0]]) #rearranging terms
A[2,0:3] = hstack([[ 0], A[2,0:2]])  
C=linalg.solve(A,B)
print ' \n The solution is \n y1 = y(0.25) = %f \n y2 = y(0.5) = %f \n y3 = y(0.75) = %f \n '%(C[0],C[1],C[2]) 

 
 The solution is 
 y1 = y(0.25) = -0.100203 
 y2 = y(0.5) = -0.137907 
 y3 = y(0.75) = -0.109079

Example No. 14_03 Pg No. 473

In [2]:
#Eigen Vectors
from numpy import eye
from sympy import symbols,Matrix,det,solve
A = [[8 ,-4],[ 2, 2 ] ]
A=Matrix(A)
lamd = symbols('lamd')
p = det(A - lamd*eye(2))
root = solve(p,lamd)
print '\n The roots are \n lamda1 = %f \n lamda2 = %f \n '%(root[0],root[1])
A1 = A - root[0]*eye(2)
X1 = Matrix([[-1*A1[0,1]/A1[0,0]],[1]])
print 'X1 = ',X1
A2 = A - root[1]*eye(2)
X2 = Matrix([[-1*A2[0,1]/A2[0,0]],[ 1]])
print 'X2 = ',X2

 The roots are 
 lamda1 = 4.000000 
 lamda2 = 6.000000 
 
X1 =  Matrix([[1.00000000000000], [1]])
X2 =  Matrix([[2.00000000000000], [1]])

Example No. 14_04 Pg No. 474

In [28]:
from numpy import array,shape,trace,poly
#Fadeev - Leverrier method

A = [[ -1, 0, 0],[ 1, -2, 3],[ 0, 2, -3 ]]
A=array(A)
r,c= shape(A)[0],shape(A)[1]
A1 = A
p= [trace(A1)]
for i in range(1,r):
    A1 = A*( A1 - p[(i-1)]*eye(3))
    p.append(trace(A1)/(i+1))
    print '\nA%d = '%(i+1)
    print A1
    print '\np%d = %f\n'%((i+1),p[i])

A2 = 
[[-5.  0.  0.]
 [ 1. -8.  9.]
 [ 0.  4. -9.]]

p2 = -11.000000


A3 = 
[[ -6.   0.   0.]
 [  1.  -6.  27.]
 [  0.   8.  -6.]]

p3 = -6.000000

Example No. 14_05 Pg No. 476

In [46]:
import numpy as np
np.seterr(divide='ignore', invalid='ignore')
from __future__ import division
from numpy.linalg import eig
from numpy import array,diag
#Eigen Vectors

A = [[ -1, 0, 0],[1, -2, 3],[0, 2, -3]]
A=array(A)
evectors,evalues = eig(A)
evectors=diag(evectors)
for i in range(0,3):
    print '\n Eigen vector - %d \n for lamda%d = %f \n X%d = '%(i+1,i+1,evalues[0,0],i+1)
    evectors[:,0] = evectors[:,0]/evectors[1,i]
    print evectors[:,i]

 Eigen vector - 1 
 for lamda1 = 0.000000 
 X1 = 
[-inf  nan  nan]

 Eigen vector - 2 
 for lamda2 = 0.000000 
 X2 = 
[ 0.  0.  0.]

 Eigen vector - 3 
 for lamda3 = 0.000000 
 X3 = 
[ 0.  0. -1.]

Example No. 14_06 Pg No. 478

In [48]:
from numpy import array, zeros,expand_dims, dot, hstack
A = array([[ 1, 2, 0],[2, 1 ,0],[0, 0, -1 ]])
X=zeros([3,8])
Y=zeros([3,7])
X[:,0] =[0,1,0]

for i in range(1,8):
    Y[:,i-1] = [xx for xx in A.dot(expand_dims(X[:,i-1], axis=1))]
    
    X[:,i] = Y[:,i-1]/max(Y[:,i-1])

print 'Iterations:\n  0      1       2      3           4            5            6             7 '
print 'X = ',X
Y=hstack([[[None],[None],[None]],Y])
print '\nY = ',Y

Iterations:
  0      1       2      3           4            5            6             7 
X =  [[ 0.          1.          0.8         1.          0.97560976  1.
   0.99726027  1.        ]
 [ 1.          0.5         1.          0.92857143  1.          0.99180328
   1.          0.99908592]
 [ 0.          0.          0.          0.          0.          0.          0.
   0.        ]]

Y =  [[None 2.0 2.0 2.8 2.8571428571428577 2.975609756097561 2.9836065573770494
  2.9972602739726026]
 [None 1.0 2.5 2.6 2.928571428571429 2.951219512195122 2.9918032786885247
  2.9945205479452053]
 [None 0.0 0.0 0.0 0.0 0.0 0.0 0.0]]