算法11.2非线性射击方法（负担和公平）Python

Question

我正在尝试编程非线性射击方法基于算法11.2从数值分析（负担和公平）。然而，运行程序后，我得到的数值结果与教科书上的答案不同。我想我的代码有问题，但我无法解决。我在图片中附上了实际的算法。Algorithm 11.2Algorithm 11.2Algorithm 11.2

这是密码

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

二阶边值问题的函数

^{pr2}$

我的结果

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

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

在fx=（1/8）*（32+2*x**3-y*yp）中，1/8的结果是0。 你应该用1/8代替。在

Answer 2

在48号线，你有

r = abs(w1[n]) - beta

而不是

^{pr2}$

做出这样的改变给出了与文本相同的解决方案

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