exponenta event banner

bvp4c

Задача решения граничных значений - метод четвертого порядка

Описание

пример

sol = bvp4c(odefun,bcfun,solinit) интегрирует систему дифференциальных уравнений вида y ′ = f (x, y), заданную odefun, в соответствии с граничными условиями, описанными bcfun и начальное предположение о решении solinit. Используйте bvpinit функция для создания начального предположения solinit, которая также определяет точки, в которых граничные условия bcfun принудительно.

пример

sol = bvp4c(odefun,bcfun,solinit,options) также использует настройки интеграции, определенные options, который является аргументом, созданным с помощью bvpset функция. Например, используйте AbsTol и RelTol для задания абсолютных и относительных допусков ошибок или FJacobian возможность предоставления аналитических частных производных odefun.

Примеры

свернуть все

Решите BVP второго порядка в MATLAB ® с помощью функций. Для этого примера используйте уравнение второго порядка

y + y = 0.

Уравнение определяется на интервале [0, δ/2] в зависимости от граничных условий

y (0) = 0,

y (δ/2) = 2.

Чтобы решить это уравнение в MATLAB, нужно записать функцию, которая представляет уравнение как систему уравнений первого порядка, функцию для граничных условий и функцию для начального предположения. Затем решатель BVP использует эти три входа для решения уравнения.

Кодовое уравнение

Напишите функцию, которая кодирует уравнение. Используйте подстановки y1 = y и y2 = y ′, чтобы переписать уравнение как систему уравнений первого порядка.

y1  = y2,

y2  =-y1.

Соответствующая функция:

function dydx = bvpfcn(x,y)
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end

Примечание.Все функции включены в конце примера как локальные функции.

Граничные условия кода

Запишите функцию, которая кодирует граничные условия в виде g (y (a), y (b)) = 0. В этой форме граничными условиями являются

y (0) = 0,

y (δ/2) -2 = 0.

Соответствующая функция:

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-2];
end

Создать начальное предположение

Используйте bvpinit чтобы создать начальное предположение для решения уравнения. Поскольку уравнение относится y к y, разумное предположение состоит в том, что решение включает тригонометрические функции. Используйте сетку из пяти точек в интервале интегрирования. Первое и последнее значения в сетке - это то, где решатель применяет граничные условия.

Функция для начального предположения принимает x в качестве входного значения и возвращает предположение для значений y1 и y2. Функция:

function g = guess(x)
g = [sin(x)
     cos(x)];
end
xmesh = linspace(0,pi/2,5);
solinit = bvpinit(xmesh, @guess);

Уравнение решения

Использовать bvp4c с функцией производной, функцией граничного условия и начальным предположением для решения задачи.

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

Решение для построения графика

plot(sol.x, sol.y, '-o')

Figure contains an axes. The axes contains 2 objects of type line.

Локальные функции

Здесь перечислены локальные функции, которые bvp4c используется для решения уравнения.

function dydx = bvpfcn(x,y) % equation to solve
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end
%--------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-2];
end
%--------------------------------
function g = guess(x) % initial guess for y and y'
g = [sin(x)
     cos(x)];
end
%--------------------------------

Решите BVP с допуском грубой ошибки с помощью двух различных решателей и сравните результаты.

Рассмотрим ОДУ второго порядка

У +2xy′+1x4y=0.

Уравнение определяется на интервале [13δ, 1] в зависимости от граничных условий

y (13δ) = 0,

y (1) = sin (1).

Чтобы решить это уравнение в MATLAB ®, необходимо записать функцию, которая представляет уравнение как систему уравнений первого порядка, записать функцию для граничных условий, задать некоторые значения опций и создать начальное предположение. Затем решатель BVP использует эти четыре входа для решения уравнения.

Кодовое уравнение

С помощью подстановок y1 = y и y2 = y ′ можно переписать ОДУ как систему уравнений первого порядка

y1  = y2,

y2′=-2xy2-1x4y1.

Соответствующая функция:

function dydx = bvpfcn(x,y)
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end

Примечание.Все функции включены в конце примера как локальные функции.

Граничные условия кода

Функция граничного условия требует, чтобы граничные условия были в форме g (y (a), y (b)) = 0. В этой форме граничными условиями являются

y (13δ) = 0,

y (1) -sin (1) = 0.

Соответствующая функция:

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-sin(1)];
end

Задать параметры

Использовать bvpset для включения отображения статистики решателя и указания допусков грубых ошибок, чтобы выделить разницу в управлении ошибками между решателями. Также для эффективности укажите аналитический якобиан

J=∂fi∂y=[∂f1∂y1∂f1∂y2∂f2∂y1∂f2∂y2]= [0 1-1x4-2x].

Соответствующей функцией, возвращающей значение якобиана, является

function dfdy = jac(x,y)
dfdy = [0      1
       -1/x^4 -2/x];
end
opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');

Создать начальное предположение

Использовать bvpinit для создания начальной догадки решения. Задайте константную функцию в качестве начального предположения с начальной сеткой из 10 точек в интервале [1/3δ, 1].

xmesh = linspace(1/(3*pi), 1, 10);
solinit = bvpinit(xmesh, [1; 1]);

Уравнение решения

Решите уравнение обоими способами bvp4c и bvp5c.

sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);
The solution was obtained on a mesh of 9 points.
The maximum residual is  9.794e-02. 
There were 157 calls to the ODE function. 
There were 28 calls to the BC function. 
sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);
The solution was obtained on a mesh of 11 points.
The maximum error is  6.742e-02. 
There were 244 calls to the ODE function. 
There were 29 calls to the BC function. 

Результаты графика

Постройте график результатов двух вычислений y1 с аналитическим решением для сравнения. Аналитическое решение:

y1 = sin (1x),

y2 =-1x2cos (1x).

xplot = linspace(1/(3*pi),1,200);
yplot = [sin(1./xplot); -cos(1./xplot)./xplot.^2];

plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo')
title('Comparison of BVP Solvers with Crude Error Tolerance')
legend('True','BVP4C','BVP5C')
xlabel('x')
ylabel('solution y')

Figure contains an axes. The axes with title Comparison of BVP Solvers with Crude Error Tolerance contains 3 objects of type line. These objects represent True, BVP4C, BVP5C.

Сюжет подтверждает, что bvp5c непосредственно управляет истинной ошибкой в вычислении, в то время как bvp4c управляет им только косвенно. При более строгих допусках ошибок эта разница между решателями не столь очевидна.

Локальные функции

Здесь перечислены локальные функции, которые решатели BVP используют для решения проблемы.

function dydx = bvpfcn(x,y) % equation to solve
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end
%---------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-sin(1)];
end
%---------------------------------
function dfdy = jac(x,y) % analytical jacobian for f
dfdy = [0      1
       -1/x^4 -2/x];
end
%---------------------------------

Входные аргументы

свернуть все

Решаемые функции, определяемые как дескриптор функции, определяющий интегрированные функции. odefun и bcfun должен принимать одинаковое количество входных аргументов.

Закодировать odefun, использовать функциональную подпись dydx = odefun(x,y) для скаляра x и вектор столбца y. Возвращаемое значение dydt - вектор столбца типа данных single или double что соответствует f (x, y ).odefun должен принимать оба входных аргумента x и y, даже если один из аргументов не используется в функции.

Например, для решения y '= 5y − 3 используйте функцию:

function dydt = odefun(t,y)
dydt = 5*y-3;
end

Для системы уравнений выход odefun является вектором. Каждый элемент вектора является решением одного уравнения. Например, для решения

y '1 = y1 + 2y2y' 2 = 3y1 + 2y2

используйте функцию:

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);
end

bvp4c также может решать задачи с сингулярностями в решении или многоточечными граничными условиями.

Пример: sol = bvp4c(@odefun, @bcfun, solinit)

Неизвестные параметры

Если решаемый BVP включает неизвестные параметры, вместо этого можно использовать функциональную подпись. dydx = odefun(x,y,p), где p - вектор значений параметров. При использовании этой функциональной сигнатуры решатель BVP вычисляет значения неизвестных параметров, начиная с начального предположения для значений параметров, представленных в solinit.

Типы данных: function_handle

Граничные условия, заданные как дескриптор функции, вычисляющий остаточную ошибку в граничных условиях. odefun и bcfun должен принимать одинаковое количество входных аргументов.

Закодировать bcfun, использовать функциональную подпись res = bcfun(ya,yb) для векторов столбцов ya и yb. Возвращаемое значение res - вектор столбца типа данных single или double соответствует остаточному значению граничных условий в граничных точках.

Например, если y (a) = 1 и y (b) = 0, то функция граничного условия

function res = bcfun(ya,yb)
res = [ya(1)-1
       yb(1)];
end
Поскольку y (a) = 1, остаточное значение ya(1)-1 должно быть 0 в точке x = a. Аналогично, поскольку y (b) = 0, остаточное значение yb(1) должно быть 0 в точке x = b.

Граничные точки x = a и x = b, где применяются граничные условия, определяются в начальной структуре приближения. solinit. Для двухточечных граничных проблем, a = solinit.x(1) и b = solinit.x(end).

Пример: sol = bvp4c(@odefun, @bcfun, solinit)

Неизвестные параметры

Если решаемый BVP включает неизвестные параметры, вместо этого можно использовать функциональную подпись. res = bcfun(ya,yb,p), где p - вектор значений параметров. При использовании этой функциональной сигнатуры решатель BVP вычисляет значения неизвестных параметров, начиная с начального предположения для значений параметров, представленных в solinit.

Типы данных: function_handle

Начальное предположение о решении, указанное как структура. Использовать bvpinit создать solinit.

В отличие от задач начального значения, задача граничного значения не может иметь решения, конечного числа решений или бесконечно большого количества решений. Важной частью процесса решения BVP является обеспечение предположения для требуемого решения. Качество этого предположения может быть критически важным для производительности решателя и даже для успешного вычисления. Некоторые рекомендации по созданию хорошего начального предположения см. в разделе Начальное угадывание решения.

Пример: sol = bvp4c(@odefun, @bcfun, solinit)

Типы данных: struct

Структура опции. Используйте bvpset для создания или изменения структуры опций.

Пример: options = bvpset('RelTol',1e-5,'Stats','on') указывает относительный допуск ошибки 1e-5 и включает отображение статистики решателя.

Типы данных: struct

Выходные аргументы

свернуть все

Структура решения. Вы можете получить доступ к полям в sol с точечным индексированием, например sol.field1. Решение (sol.x,sol.y) является непрерывным на интервале интегрирования, определенном в начальной сетке solinit.x и имеет непрерывную первую производную. Вы можете использовать sol с deval для оценки решения в других точках интервала.

Структура sol имеет эти поля.

ОбластьОписание

x

Сетка, выбранная bvp4c. Эта сетка обычно содержит точки, отличные от начальной сетки solinit.x.

y

Приближение к y (x) в точках сетки sol.x.

yp

Приближение к y(x) в точках сетки sol.x.

parameters

Окончательные значения неизвестных параметров, указанных в solinit.parameters.

solver

'bvp4c'

stats

Статистика вычислительных затрат, связанная с решением: количество точек ячеистой сети, остаточная ошибка и количество вызовов odefun и bcfun.

Подробнее

свернуть все

Проблемы с многоточечными граничными значениями

Для проблем с многоточечными граничными значениями граничные условия применяются в нескольких точках в интервале интегрирования.

bvp4c может решать задачи многоточечных граничных значений, где a = a0  < a1  < a2  < ... < an  = b - граничные точки в интервале [a, b]. Точки a1, a2,..., a − 1 представляют интерфейсы, которые делят [a, b] на области.bvp4c перечисляет области слева направо (от a до b) с индексами, начинающимися с 1. В регионе k, [ak 1, ak],bvp4c оценивает производную как

yp = odefun(x,y,k) 

В функции граничных условий bcfun(yleft,yright), yleft(:,k) - решение на левой границе [ak 1, ak]. Аналогично,yright(:,k) - решение на правой границе области k. В частности, yleft(:,1) = y(a) и yright(:,end) = y(b).

При создании начального предположения с помощью bvpinit, используйте двойные записи в xinit для каждой точки интерфейса. См. справочную страницу для bvpinit для получения дополнительной информации.

Если yinit является функцией, bvpinit требования y = yinit(x,k) чтобы получить начальное предположение для решения на x в регионе k. В структуре решения sol возвращено bpv4c, sol.x имеет двойные записи для каждой точки интерфейса. Соответствующие столбцы sol.y содержат левое и правое решения на интерфейсе соответственно.

Пример решения проблемы с трехточечным граничным значением см. в разделе Решение BVP с несколькими граничными условиями.

Проблемы с сингулярными граничными значениями

bvp4c решает класс проблем с сингулярными граничными значениями, включая проблемы с неизвестными параметрами p, формы

y '= Syx + f (x, y, p), 0 = bc (y (0), y (b), p).

Интервал должен быть [0, b] с b > 0. Часто такие проблемы возникают при вычислении гладкого решения ОДУ, которые являются результатом дифференциальных уравнений в частных производных (ПДЭ) из-за цилиндрической или сферической симметрии. Для сингулярных задач задается матрица (константа). S как значение 'SingularTerm' вариант bvpset, и odefun вычисляет только f (x, y, p). Граничные условия и начальное предположение должны соответствовать необходимому условию для гладкости S· y (0) = 0.

Пример решения проблемы сингулярного граничного значения см. в разделе Решение BVP с сингулярным термином.

Алгоритмы

bvp4c - код конечных разностей, реализующий трёхступенчатую формулу Лобатто Λ a [1], [2]. Это формула словосочетания, и многочлен словосочетания обеспечивает C1-continuous решение, которое является точным четвёртого порядка равномерно в интервале интегрирования. Выбор сетки и управление ошибками основаны на остатке непрерывного решения.

Метод коллокации использует сетку точек для разделения интервала интегрирования на субинтервалы. Решатель определяет числовое решение, решая глобальную систему алгебраических уравнений, вытекающих из граничных условий, и условий коллокации, наложенных на все субинтервалы. Затем решатель оценивает ошибку числового решения на каждом субинтервале. Если решение не удовлетворяет критериям допуска, решатель адаптирует сетку и повторяет процесс. Необходимо указать точки начальной сетки, а также исходную аппроксимацию решения в точках сетки.

Ссылки

[1] Шампин, Л.Ф. и Дж. Кьерзенка. «Решатель BVP, основанный на остаточном контроле и MATLAB PSE». ACM Trans. Math. Softw. Том 27, номер 3, 2001, стр. 299-316.

[2] Шампин, Л.Ф., М.У. Райхельт и Ж. Кьерзенка. «Решение граничных задач для обычных дифференциальных уравнений в MATLAB с bvp4c». Обмен файлами MATLAB, 2004.

Представлен до R2006a