Algorithm 11.2 Nonlinear Shooting Method (Burden and Faires) Python

Question

I am trying to program nonlinear shooting method based on algorithm 11.2 from Numerical Analysis (Burden and Faires). However, after running the program, the numerical result i am getting is different from the answer in the textbook. i think there is something wrong in my coding but i cannot figure it out. I attached actual algorithm in the picture. Algorithm 11.2 Algorithm 11.2 Algorithm 11.2

Here is the code

from numpy import zeros, abs

def shoot_nonlinear(a,b,alpha, beta, n, tol, M):

    w1 = zeros(n+1)  
    w2 = zeros(n+1)
    h = (b-a)/n
    k = 1
    TK = (beta - alpha)/(b - a)

    print("i""  x" "     " "W1""     " "W2")
    while k <= M:

        w1[0] = alpha
        w2[0] = TK
        u1    = 0
        u2    = 1

        for i in range(1,n+1):
            x = a + (i-1)*h     #step 5

            t = x + 0.5*(h)

            k11 = h*w2[i-1]     #step 6

            k12 = h*f(x,w1[i-1],w2[i-1])
            k21 = h*(w2[i-1] + (1/2)*k12)
            k22 = h*f(t, w1[i-1] + (1/2)*k11, w2[i-1] + (1/2)*k12)
            k31 = h*(w2[i-1] + (1/2)*k22)
            k32 = h*f(t, w1[i-1] + (1/2)*k21, w2[i-1] + (1/2)*k22)
            t   = x + h
            k41 = h*(w2[i-1]+k32)
            k42 = h*f(t, w1[i-1] + k31, w2[i-1] + k32)
            w1[i] = w1[i-1] + (k11 + 2*k21 + 2*k31 + k41)/6
            w2[i] = w2[i-1] + (k12 + 2*k22 + 2*k32 + k42)/6   
            kp11 = h*u2
            kp12 = h*(fy(x,w1[i-1],w2[i-1])*u1 + fyp(x,w1[i-1], w2[i-1])*u2)
            t    = x + 0.5*(h)
            kp21 = h*(u2 + (1/2)*kp12)
            kp22 = h*((fy(t, w1[i-1],w2[i-1])*(u1 + (1/2)*kp11)) + fyp(x+h/2, w1[i-1],w2[i-1])*(u2 +(1/2)*kp12))
            kp31 = h*(u2 + (1/2)*kp22)
            kp32 = h*((fy(t, w1[i-1],w2[i-1])*(u1 + (1/2)*kp21)) + fyp(x+h/2, w1[i-1],w2[i-1])*(u2 +(1/2)*kp22))
            t    = x + h
            kp41 = h*(u2 + kp32)
            kp42 = h*(fy(t, w1[i-1], w2[i-1])*(u1+kp31) + fyp(x + h, w1[i-1], w2[i-1])*(u2 + kp32))
            u1 = u1 + (1/6)*(kp11 + 2*kp21 + 2*kp31 + kp41)
            u2 = u2 + (1/6)*(kp12 + 2*kp22 + 2*kp32 + kp42)


        r = abs(w1[n]) - beta
        #print(r)
        if r < tol:
            for i in range(0,n+1):
                x = a + i*h
                print("%.2f %.2f %.4f %.4f" %(i,x,w1[i],w2[i]))
            return

        TK = TK -(w1[n]-beta)/u1

        k = k+1


    print("Maximum number of iterations exceeded")   
    return

function for 2nd order boundary value problem

def f(x,y,yp):
fx = (1/8)*(32 + 2*x**3 -y*yp)
return fx

def fy(xp,z,zp):
fyy = -(1/8)*(zp)
return fyy

def fyp(xpp,zpp,zppp):
fypp = -(1/8)*(zpp)
return fypp

a = 1         # start point
b = 3         # end point
alpha = 17    # boundary condition
beta = 43/3   # boundary condition
N = 20        # number of subintervals
M = 10        # maximum number of iterations
tol = 0.00001 # tolerance


shoot_nonlinear(a,b,alpha,beta,N,tol,M)

My result

i      x     W1      W2
0.00 1.00 17.0000 -16.2058
1.00 1.10 15.5557 -12.8379
2.00 1.20 14.4067 -10.2482
3.00 1.30 13.4882 -8.1979
4.00 1.40 12.7544 -6.5327
5.00 1.50 12.1723 -5.1496
6.00 1.60 11.7175 -3.9773
7.00 1.70 11.3715 -2.9656
8.00 1.80 11.1203 -2.0783
9.00 1.90 10.9526 -1.2886
10.00 2.00 10.8600 -0.5768
11.00 2.10 10.8352 0.0723
12.00 2.20 10.8727 0.6700
13.00 2.30 10.9678 1.2251
14.00 2.40 11.1165 1.7444
15.00 2.50 11.3157 2.2331
16.00 2.60 11.5623 2.6951
17.00 2.70 11.8539 3.1337
18.00 2.80 12.1883 3.5513
19.00 2.90 12.5635 3.9498
20.00 3.00 12.9777 4.3306

Actual result for w1

x     W1
1.0   17.0000
1.1   15.7555
1.2   14.7734
1.3   13.3886
1.4   12.9167
1.5   12.5601
1.6   12.3018
1.7   12.1289
1.8   12.0311
1.9   12.0000
2.0   12.0291
2.1   12.1127
2.2   12.2465
2.3   12.4267 
2.4   12.6500
2.5   12.9139
2.6   13.2159
2.7   13.5543
2.8   13.9272
2.9   14.3333
3.0   14.7713

Answer 1

On Line 48, you have

r = abs(w1[n]) - beta

Instead of

r = abs(w1[n] - beta)

Making this change gives the same solution as the text,

x     W1
1.0   17.0000
1.1   15.7555
1.2   14.7734
1.3   13.9978
1.4   13.3886
1.5   12.9167
1.6   12.5601
1.7   12.3018
1.8   12.1289
1.9   12.0311
2.0   12.0000
2.1   12.0291
2.2   12.1127
2.3   12.2465
2.4   12.4267 
2.5   12.6500
2.6   12.9139
2.7   13.2159
2.8   13.5543
2.9   13.9272
3.0   14.3333

Answer 2

As a remark in fx = (1/8)*(32 + 2*x**3 -y*yp), 1/8 is going to give the result 0. You should use 1./8 instead.