Для решения краевых задач, т.е. дифференциальных уравнений произвольного порядка n, или систем из n уравнений первого порядка вида , служит РЕШАТЕЛЬ bvp4c.
Точность решения контролируется при помощи двух опций управляющей структуры RelTol – для относительной точности и AbsTol – для абсолютной. Сама управляющая структура генерируется функцией bvpset
и задается в качестве четвертого входного аргумента солвера bvp4c.
Рассмотрим решение краевых задач на примере обыкновенного дифференциального уравнения второго порядка (n=2). Требуется найти функцию у(х), удовлетворяющую на отрезке [a, b] дифференциальному уравнению и граничными условиями:
, (первое граничное условие в точке a)
, (второе граничное условие в точке b)
где – неизвестные функции, а – заданные числа.
Решение этой задачи состоит из следующих этапов:
1. Преобразование дифференциального уравнения второго порядка к системе двух уравнений первого порядка.
2. Написание функции для вычисления правой части системы.
4. Формирование начального приближения при помощи специальной функции bvpinit.
5. Вызов солвера bvp4c для решения граничной задачи.
6. Визуализация результата.
Этап 1. Введение вспомогательных функций y1(x)=y и y2(x)=y1¢=y¢, приводит к системе из двух уравнений первого порядка относительно них (при постановке граничных условий для вспомогательных функций требуется преобразовать их так, чтобы в правых частях стояли нули):
,
,
, (первое граничное условие в точке a)
(второе граничное условие в точке b).
Этап 2. Функция правой части системы зависит от х и вектора у, состоящего из двух компонент, у(1) соответствует , а у(2) – .
Этап 3. Функция, описывающая граничные условия, зависит от двух аргументов векторов ya и yb и имеет следующую структуру:
function g = bound(ya,yb)
g = [alpha1*ya(1)+beta1*ya(2)–gamma1; alpha2*yb(1)+beta2*yb(2) –gamma2];
причем вместо alpha1, beta1, gamma1, alpha2, beta2, gamma2 следует подставить заданные числа.
Этап 4. Солвер bvp4c основан на методе конечных разностей, т. е. получающееся решение есть векторы значений неизвестных функций в точках отрезка (в узлах сетки). Выбор начальной сетки и приближение может оказать влияние на решение, выдаваемое солвером bvp4c. Для задания начальной сетки и приближения предназначена функция bvpinit, обращение к которой в самом простом случае выглядит в виде
initsol = bvpinit (начальная сетка, начальное приближение к решению),
где начальная сетка определяется вектором координат узлов, упорядоченных по возрастанию и принадлежащих отрезку [a, b]. Если имеется априорная информация о решении, то разумно среди точек начальной сетки указать те, в которых решение сильно изменяется.
Формирование равномерной сетки целесообразно производить функцией linspace,которая возвращает вектор х из n равноотстоящих узлов между a и b, включая границы:
X = linspace (a, b, n).
Заданная сетка модифицируется солвером в процессе решения для обеспечения требуемой точности. Постоянное начальное приближение задается вектором из двух элементов функций . Начальное приближение может зависеть от х, в этом случае:
1) требуется запрограммировать функцию, которая по заданному х возвращает вектор из двух компонент со значениями ,
2) поместить указатель на нее во втором входном аргументе bvpinit. В результате работы bvpinit генерируется структура initsol с информацией о начальном приближении.
Этап 5. После определения начального приближения вызывается солвер bvp4c, входными аргументами которого являются указатели на функции правой части системы и граничных условий, начальное приближение и, при необходимости, управляющая структура с параметрами вычислительного процесса. Управляющая структура формируется при помощи функции bvpset. Солвер bvp4c возвращает единственный выходной аргумент – структуру с информацией о расчетной сетке, значения искомой функции и ее производной. Солвер bvp4c так же, как и и солвер dde23 для ДУ с запаздыванием, позволяет получить результат только в виде структуры.
Для вычисления значений приближенного решения в произвольных точках отрезка следует применить функцию deval, так же как и в случае решения задачи Коши.
Рассмотрим решение граничной задачи при помощи bvp4c для дифференциального уравнения второго порядка:
при граничных условиях
После преобразования получаем следующую систему ДУ первого порядка (этап 1)
,
,
,
.
Подготовим функцию по имени rside для правых частей системы ДУ в виде m-файла (этап 2).
function f = rside(x, y)
f = [y(2); –sin(x)];
Подготовим функцию по имени bound дл я граничных условий в виде m-файла (этап 3).
function f = bound(ya, yb)
f = [ya(1); yb(1) + yb(2) + 1];
При помощи bvpinit зададим начальную сетку на отрезке [0, ], например из 20 узлов, и постоянное начальное приближение: , .
% формирование начальной сетки meshinit и начального приближения yinit
»meshinit = linspace(0, 11*pi/2, 20);
»yinit = [1 0];
»initsol = bvpinit(meshinit, yinit);
Вызовем солвер bvp4c и получим приближенное решение.
»sol = bvp4c(@rside, @bound, initsol);
Отобразим графически приближенное решение, извлекая нужные компоненты из структуры возвращаемой солвером, т.е. используем поля x и y структуры sol для построения решения:
% в sol.x содержатся координаты сетки
% в sol.y матрица
% в sol.y(1, :) соответствует значениям функции y1 в sol.x(:)
% в sol.y(2, :) соответствует значениям функции y2 в sol.x(:)
»plot(sol.x, sol.y(1, :), 'k.')
Отобразим графически точное решение на отрезке [0, ] для сравнения.
»hold on
»fplot(@sin, [0 11*pi/2]);
»title('Решение граничной задачи при помощи bvp4c');
Как следует из полученных графиков (в виде точек и сплошной линии) для простых задач солвер bvp4c дает хорошие результаты. По умолчанию решение найдено с относительной точностью 10-3 и относительной 10-6 по невязке.
, 
служит РЕШАТЕЛЬ bvp4c.
, (первое граничное условие в точке a)
, (второе граничное условие в точке b)
– неизвестные функции, а 
– заданные числа.
,
,
, (первое граничное условие в точке a)
(второе граничное условие в точке b).
, а у(2) – 
.
. Начальное приближение может зависеть от х, в этом случае:
,
при граничных условиях 
,
,
.
], например из 20 узлов, и постоянное начальное приближение: 
, 
.
на отрезке [0, 