Метод стрельбы решения краевой задачи

import math
from array import array
 
def y1(x):
    return math.e ** x + math.e ** (2 * x) + (0.1 * x - 0.12) * math.cos(x) - (0.3 * x + 0.34) * math.sin(x)
 
def g(x,y,z):
    return x * math.cos(x) + 3 * z - 2 * y
    
a,b = 0.0,1.0 
n = 10.0
eps = 0.001
 
def W(ksi):
    x=array('f')
    y=array('f')
    z=array('f')
 
    x.insert(0,a)
    y.insert(0,ksi)
    z.insert(0,2.76)
    i = 0
    h = (b - a) / n
    while i < int(n-1):
        x.append(x[i] + h)
        y.append(y[i] + h * z[i])
        z.append(z[i] + h * g(x[i],y[i],z[i]))
        i += 1    
    return y[int(n-1)]
 
 
min_ksi = -1000000.0
max_ksi = 1000000.0
 
while True:
    ksi = (min_ksi + max_ksi) / 2
    my_result = W(ksi)
    goal = y1(b)
    if math.fabs(my_result-goal) < eps:
        break
    if my_result > goal:
        min_ksi = ksi
    else:
        max_ksi = ksi
 
print(ksi, ' ',W(ksi),' ',y1(b))

def newton_method(f, df, x0, eps):
    """Ищет решение f(x)=0 методом Ньютона
    
    f - функция от x, для которой ищем решение
    df - производная по x от f
    x0 - начальная точка
    eps - погрешность
    
    returns:
    x, для которого abs(f(x)) < eps
    """
    
    while True:
        delta = abs(f(x0))
        if delta < eps:
            return x0
        x0 = x0 - f(x0)/df(x0)
 
 
def shoot(fun_y2, y0, ksi, x0=0.0, x1=1.0, n=1000):
    """Рассчитывает один "выстрел" по методу стрельбы
    
    fun_y2 - вторая производная функции y
    y0 - y(x0)
    ksi - "угол выстрела", т.е. наклон касательной к y, т.е. y(x0)
    x0, x1 - диапазон x, на котором производится расчет
    n - количество блоков, на котороые разбивается диапазон (x0, x1) для расчета
    
    returns:
    конечное значение пули после выстрела, т.е. y(x1)
    
    """
    x = x0
    y = y0
    y1 = ksi
    
    delta = (x1 - x0) / n
    for i in range(n):
        x += delta
        y += delta * y1
        y1 += delta * fun_y2(x,y,y1)
    return y
 
 
def shooting_method(fun_ode, x0, y0, x1, y1, eps=1e-3):
    """Поиск решение ОДУ методом стрельбы
    
    fun_ode(x,y,y1) - p(x), в выражении если y''=p(x,y,y')
    (x0,y0) - начальная точка
    (x1,y1) - конечная точка
    1e-3 - погрешность
    
    return:
    угол стрельбы, y'(x0)
    """
    
    def F(ksi):
        return shoot(fun_ode, y0, ksi, x0, x1) - y1
    
    def dF(ksi):
        return F(ksi+eps) - F(ksi) / eps
    
    return newton_method(F, dF, x0, eps)

# Проверка метода Ньютона. Решим уравнение x*x-100 = 0
newton_method(lambda x: x*x-100, lambda x: 2*x, 5, eps=1e-6)
# 10.00000000000001
 
 
 
# Проверка выстрела
 
def fun_y2(x, y, y1):
    return x**4
 
def fun_y1(x, c1):
    return x**5 / 5 + c1
 
# Для любых c1 и c2 не должен срабатывать assert
c1=1
c2=1
print("c1=", c1)
print("c2=", c2)
 
x0 = 0.0
x1 = 1.0
y0 = fun_y(x0, c1=c1, c2=c2)
y1 = fun_y(x1, c1=c1, c2=c2)
ksi = fun_y1(x0, c1)
 
calculated_y1 = shoot(fun_y2, y0, ksi)
print("y1=", y1)
print("расчетный y1=", calculated_y1)
assert abs(calculated_y1 - y1) < 1e-3, "разница должна быть мала для любых c1 и c2"
#c1= 1
#c2= 1
#y1= 2.033333333333333
#расчетный y1= 2.0333332500000516
 
 
 
# Проверка краевой задачи
def fun_ode(x, y, y1):
    return x**4
 
def fun_y1(x, c1):
    return x**5 / 5 + c1
 
def fun_y(x, c1, c2):
    return x**6 / 30 + c1 * x + c2
 
# Для любых c1 и c2 не должен срабатывать assert
c1=1
c2=1
print("c1=", c1)
print("c2=", c2)
 
x0 = 0.0
x1 = 1.0
y0 = fun_y(x0, c1=c1, c2=c2)
y1 = fun_y(x1, c1=c1, c2=c2)
ksi = fun_y1(x0, c1)
 
calculated_ksi = shooting_method(fun_ode, x0, y0, x1, y1)
print("ksi=", ksi)
print("расчетный ksi=", calculated_ksi)
assert abs(calculated_ksi - ksi) < 1e-3, "разница должна быть мала для любых c1 и c2"
 
#c1= 1
#c2= 1
#ksi= 1.0
#расчетный ksi= 0.9999925196666988