Метод стрельбы

Метод стрельбы (краевая задача) — численный метод, заключающийся в сведении краевой задачи к некоторой задаче Коши для той же системы дифференциальных уравнений. Суть: первое решение при последовательном изменении аргумента и повторении вычислений становится точнее

Описание методаПравить

Рассматривается задача для системы двух уравнений первого порядка с краевыми условиями общего вида:

система

$u^{\prime }(x)=f(x,u,v)$$u'(x)=f(x,u,v)$ 
$v^{\prime }(x)=g(x,u,v)$$v'(x)=g(x,u,v)$ 

граничные условия

$a⩽<!--\; ⩽\; -->x⩽<!--\; ⩽\; -->b$$a\leqslant x\leqslant b$ 
$\mathrm{\phi <!--\; \phi \; -->}[u(a),v(a)]=0$$\varphi [u(a),v(a)]=0$ 
$\mathrm{\psi <!--\; \psi \; -->}[u(b),v(b)]=0$$\psi [u(b),v(b)]=0$ 

АлгоритмПравить

1. Выбирается произвольно условие $u(a)=\mathrm{\eta <!--\; \eta \; -->}$$u(a)=\eta$ .

2. Рассматривается левое краевое условие как алгебраическое уравнение $\mathrm{\phi <!--\; \phi \; -->}(\mathrm{\eta <!--\; \eta \; -->},v(a))=0$$\varphi (\eta ,v(a))=0$ . Определяем удовлетворяющее ему значение $v(a)=\mathrm{\zeta <!--\; \zeta \; -->}(\mathrm{\eta <!--\; \eta \; -->})$$v(a)=\zeta (\eta )$ .

3. Выбираются значения $u(a)=\mathrm{\eta <!--\; \eta \; -->},v(a)=\mathrm{\zeta <!--\; \zeta \; -->}$$u(a)=\eta ,v(a)=\zeta$  в качестве начальных условий задачи Коши для рассматриваемой системы и интегрируется эта задача Коши любым численным методом (например, по схемам Рунге — Кутты).

4. В итоге получается решение $u(x;\mathrm{\eta <!--\; \eta \; -->}),v(x;\mathrm{\eta <!--\; \eta \; -->})$$u(x;\eta ),v(x;\eta )$ , зависящее от η как от параметра.

Значение $\mathrm{\zeta <!--\; \zeta \; -->}$$\zeta$  выбрано так, что найденное решение удовлетворяет левому краевому условию. Однако правому краевому условию это решение, вообще говоря, не удовлетворяет: при его подстановке левая часть правого краевого условия, рассматриваемая как некоторая функция параметра $\mathrm{\eta <!--\; \eta \; -->}$$\eta$ :

$\stackrel{~<!--\; ~\; -->}{\mathrm{\psi <!--\; \psi \; -->}}(\mathrm{\eta <!--\; \eta \; -->})=\mathrm{\psi <!--\; \psi \; -->}(u(b;\mathrm{\eta <!--\; \eta \; -->}),v(b;\mathrm{\eta <!--\; \eta \; -->}))$${\tilde {\psi }}(\eta )=\psi (u(b;\eta ),v(b;\eta ))$ ,

не обратится в нуль.

5. Подбирается параметр η по условию нахождения такого значения, для которого $\stackrel{~<!--\; ~\; -->}{\mathrm{\psi <!--\; \psi \; -->}}(\mathrm{\eta <!--\; \eta \; -->})\approx <!--\; \approx \; -->0$${\tilde {\psi }}(\eta )\approx 0$  с требуемой точностью.

Таким образом, решение краевой задачи сводится к нахождению корня одного алгебраического уравнения $\stackrel{~<!--\; ~\; -->}{\mathrm{\psi <!--\; \psi \; -->}}(\mathrm{\eta <!--\; \eta \; -->})=0$${\tilde {\psi }}(\eta )=0$ .

Пример программы на языке PythonПравить

import matplotlib.pyplot as plt
import numpy as np

a, b = 0.0, 1.0
A, B = 1.0, np.e
n = 5
h = (b - a) / n
D0, D1 = A + h, h

y = [[A, D0], [0, D1]]

def p(x):   return 1

def q(x):   return 1

def f(x):   return 3 * (np.e **x)

def get_c1():
    global n
    return (B - y[0][n]) / y[1][n]

def get_solv_y_i(i): return y[0][i] + get_c1() * y[1][i]

x = np.linspace(a, b, n+1)

def div(a, b):
    return a / b

for i in range(1, n):
    y[0].append(
        div(
            (h ** 2 * f(x[i]) - (1.0 - (h / 2) * p(x[i])) * y[0][i - 1] - (h ** 2 * q(x[i]) - 2) * y[0][i]),
            1 + h / 2 * p(x[i])
        )
    )
    y[1].append(
        div(
            -(1 - h / 2 * p(x[i])) * y[1][i - 1] - (h ** 2 * q(x[i]) - 2) * y[1][i],
            1 + h / 2 * p(x[i])
        )
    )

plt.plot(x, [get_solv_y_i(i) for i in range(n + 1)])
plt.show()

for i in range(n):
    print(x[i], get_solv_y_i(i))