Я пытаюсь запрограммировать метод нелинейной съемки, основанный на алгоритме 11.2 из численного анализа (бремя и Faires). Однако после запуска программы числовой результат, который я получаю, отличается от ответа в учебнике. я думаю, что что-то не так в моем кодировании, но я не могу понять это. Я приложил фактический алгоритм на картинке. Алгоритм 11.2
Алгоритм 11.2
Алгоритм 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
функция для краевой задачи 2-го порядка
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