问算法11.2非线性射击方法(负担和计算) Python

EN

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

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)

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

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