3.3 Решение краевой задачи методом стрельбы

Рассмотрим краевую задачу для уравнения второго порядка, разрешенного относительно второй производной:

y = f(x, y, y). (3.3)

Будем искать решение y = y(x) этого уравнения на отрезке [a,b]. Граничные условия на концах рассматриваемого отрезка примем в простом виде

y(a)=y₀, y(b)=y₁, (3.4)

что получается из (3.2) при ₀=₁=0 и соответствует граничным условиям 1-го рода.

Сущность метода стрельбы заключается в сведении решения краевой задачи (3.3), (3.4) к решению последовательности задач Коши для того же уравнения (3.3) с начальными условиями

y(a) = y₀, y(a) = tg . (3.5)

Здесь y₀ – точка на оси ординат, в которой помещается начало искомой интегральной кривой;  – угол наклона касательной к интегральной кривой в этой точке. Считая решение задачи Коши y=y(x, ) зависящим от параметра , будем искать такую интегральную кривую y=y(x, *), которая выходит из точки (a, y0) и попадает в точку (b, y₁). Таким образом, если  = *, то решение y(x, ) задачи Коши совпадает с решением y(x) краевой задачи. При х = b, учитывая второе граничное условие (3.4), получаем y(b, ) = y₁, или

y(b, ) – y₁=0. (3.6)

Следовательно, для нахождения параметра  получим уравнение вида F() = 0, где F() = y(b, ) – y₁. Это уравнение отличается от привычной записи тем, что функцию F() нельзя представить в виде некоторого аналитического выражения, поскольку она является решением задачи Коши (3.3), (3.4). Тем не менее, для решения уравнения (3.6) может быть использован любой из рассмотренных ранее методов решения нелинейных уравнений. В частности, одним из самых надежных является метод Ньютона. Его применение состоит в следующем. Пусть ₀ – начальное приближение к *. Построим итерационный процесс для нахождения последующих приближений _k с помощью формулы метода Ньютона:

.

С учетом того, что , имеем

. (3.7)

Производную в знаменателе этого выражения можно найти численно:

. (3.8)

Здесь  – произвольное малое возмущение .

Для вычисления правой части (3.8) нужно решить задачу Коши при  = _k–1+, в результате чего найдем значение y(x, _k–1+). Затем по формуле (3.7) находим следующее приближение _k параметра * и т.д. Этот итерационный процесс продолжается до тех пор, пока два последовательных приближения _k–1 и _k не станут отличаться меньше, чем на заданное малое число .

Описанный алгоритм называется методом стрельбы вполне оправданно, поскольку в нем как бы проводится «пристрелка» по углу наклона интегральной кривой в начальной точке с тем, чтобы попасть в точку (b, y₁). Следует отметить, что этот алгоритм хорошо работает в том случае, если решение y(x, ) не слишком чувствительно к изменениям ; в противном случае можно столкнуться с неустойчивостью.

Замечание. Вместо угла  итерационный процесс можно строить для тангенса . Вычислительная процедура при этом почти не меняется.

Пример 4. Решить краевую задачу:

, , .

Решение. Имеем краевую задачу 1-го рода. Программа на языке Maple может быть реализована следующим образом.

> restart; # очистка оперативной памяти

> eqn:=(D@@2)(y)(x)-y(x)=x; # задание дифференциального уравнения

> tgphi:=0: eps:=0.001: # инициализация переменных

> n:=0; Y1:=0; YB:=0; a:=0; b:=1;

Поясним смысл введенных переменных: tgphi – tg – тангенс угла наклона кривой в левой граничной точке x=0, именно эта величина подлежит уточнению с помощью итераций; eps – малое число – используется в условии прекращения итераций; а также для приращения величины tg; n – число итераций, Y1 – точное значение y₁ из граничного условия в правой точке: y(1)=y₁; YB – приближение y₁; a, b – левая и правая граничные точки.

> Digits:=15; # задание числа значащих цифр в вычислениях (по умолчанию 10)

Далее идет цикл, на каждом шаге которого решаются две близкие задачи Коши, отличающиеся только граничным условием в точке x=1, и вычисляется поправка к tgphi. Как только значение текущей задачи Коши при x=1 (переменная YB) приблизится к y₁, цикл останавливается.

> while (abs(YB-Y1)>eps) or (n=0) do

> p0:=dsolve({eqn,y(a)=0.,D(y)(a)=tgphi},y(x),

type=numeric,output=listprocedure); # решение началь­ной задачи с условиями y(0)=0, y(0)=tg

> f:=subs(p0,y(x));YB:=f(b); # YB – значение найденного решения при x=1

> p1:=dsolve({eqn,y(a)=0.,D(y)(a)=tgphi+eps},

y(x),type=numeric,output=listprocedure); # решение возмущенной начальной задачи с условиями y(0)=0, y(0)=tg+

> f1:=subs(p1,y(x)); # YB – значение возмущенного решения при x=1

> n:=n+1; # счетчик числа итераций

> tgphi:=tgphi+(Y1-YB)/(f1(b)-YB)*eps; # уточнение tg по методу Ньютона

> od;

> tgphi; YB; # вывод tg, равного y(0), и YB – значения y(x) в правой граничной точке x=1

> g1:=plots[odeplot](p0,a..b,colour=black):

> plots[display](g1);

Интересно сравнить полученный результат с точным решением, тем более, что данная краевая задача имеет аналитическое решение, которое без труда найдем с помощью встроенной функции dsolve:

> p:=dsolve({eqn,y(0)=0.,y(1)=0},y(x));

> u:=unapply(subs(p,y(x)),x);

Протабулируем полученное методом стрельбы решение вместе с точным решением с шагом 0.1:

> for t from a to b by 0.1 do

> printf(`x=%g y=%g u=%g\n`,t,f(t),evalf(u(t)));

> od;

> g2:=plot(u(x),x=a..b): # график точного решения

> plots[display](g1,g2); # совмещение графиков приближен­ного и точного решений

Решение граничной задачи второго рода

y = f(x, y, y), y(a) = M₀, y(b) = M₁

по методу стрельбы сводится к последовательному решению задач Коши

y = f(x, y, y), y(a) = M₀, y(a) = y₀,

причем «пристрелка» проводится по параметру y₀, значение которого целенаправленно подбирается так, чтобы найденное решение приближенно удовлетворяло условию в точке x=b: .

Пример 5. Решить краевую задачу

, , .

Решение в системе компьютерной математики Maple:

> restart;

> Digits:=15; eps:=0.001; a:=0; b:=1; M0:=-1.; M1:=1; # eps – малое число – используется в условии завершения цикла, а также для приращения величины y₀; a, b – левая и правая граничные точки; M0 и M1 – заданные граничные значения y на концах отрезка [a, b].

> n:=0; y0:=0; # n – число итераций, y0 – приближение величины y₀

> eqn:=(D@@2)(y)(x)-y(x)=x;

> while (abs(dy1-M1)>eps) or (n=0) do

> p0:=dsolve({eqn,y(a)=y0,D(y)(a)=M0},y(x),

type=numeric,output=listprocedure);

> q:=subs(p0,diff(y(x),x)); dy1:=q(b);

> p1:=dsolve({eqn,y(a)=y0+eps,D(y)(a)=M0},y(x),

type=numeric,output=listprocedure);

> q1:=subs(p1,diff(y(x),x));

> n:=n+1;

> y0:=y0-(dy1-M1)/(q1(b)-dy1)*eps;

> od;

> y0; dy1; # dy1 – приближение величины M₁

> g1:=plots[odeplot](p0,a..b,colour=black):

> plots[display](g1);

Рис. 3.2.

> for t from a to b by 0.1 do printf(`x=%g

y=%g \n`,t,subs((p0[2]),y(x))(t)); od;