算法 11.2 非线性射击方法（负担和费尔斯）Python
Algorithm 11.2 Nonlinear Shooting Method (Burden and Faires) Python
我正在尝试根据数值分析（Burden 和 Faires）中的算法 11.2 对非线性射击方法进行编程。然而，在运行程序之后，我得到的数值结果与教科书上的答案不同。我认为我的编码有问题，但我无法弄清楚。我在图片中附上了实际的算法。 Algorithm 11.2
Algorithm 11.2
Algorithm 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
二阶边值问题的函数
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
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
如 fx = (1/8)*(32 + 2*x**3 -y*yp) 中的注释，1/8 将给出结果 0。
您应该改用 1./8。
在第 48 行，您有
r = abs(w1[n]) - beta
而不是
r = abs(w1[n] - beta)
进行此更改会得到与文本相同的解决方案，
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
我正在尝试根据数值分析（Burden 和 Faires）中的算法 11.2 对非线性射击方法进行编程。然而，在运行程序之后，我得到的数值结果与教科书上的答案不同。我认为我的编码有问题，但我无法弄清楚。我在图片中附上了实际的算法。 Algorithm 11.2 Algorithm 11.2 Algorithm 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
二阶边值问题的函数
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
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
如 fx = (1/8)*(32 + 2*x**3 -y*yp) 中的注释，1/8 将给出结果 0。 您应该改用 1./8。
在第 48 行，您有
r = abs(w1[n]) - beta
而不是
r = abs(w1[n] - beta)
进行此更改会得到与文本相同的解决方案，
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