# Chapter 14 - Boundary value and eigen value problems¶

## Example No. 14_01 Pg No. 467¶

In :
#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 = z[i]
m1 = f(x[i],y[i])
m2 = z[i] + h*m1
m2 = f(x[i]+h,y[i]+h*m1)
m = (m1 + m2)/2
m = ( m1 + m2 )/2
x.append(x[(i)] + h)
y.append(y[(i)] + h*m)
z.append(z[(i)] + h*m)

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 :
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], ]) #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,C,C)


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 :
#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,root)
A1 = A - root*eye(2)
X1 = Matrix([[-1*A1[0,1]/A1[0,0]],])
print 'X1 = ',X1
A2 = A - root*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], ])
X2 =  Matrix([[2.00000000000000], ])


## Example No. 14_04 Pg No. 474¶

In :
from numpy import array,shape,trace,poly

A = [[ -1, 0, 0],[ 1, -2, 3],[ 0, 2, -3 ]]
A=array(A)
r,c= shape(A),shape(A)
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 :
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 :
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]]