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, разумное предположение - то, что решение включает тригонометрические функции. Используйте сетку пяти точек в интервале интегрирования. Первые и последние значения в mesh - то, где решатель применяет граничные условия.

Функция для исходного предположения принимает 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')

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

Перечисленный здесь локальные функции что 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 в допуске грубой ошибки с двумя другими решателями и сравните результаты.

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

y+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=fiy=[f1y1f1y2f2y1f2y2]=[01-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')

График подтверждает тот 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'=5y3, используйте функцию:

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. Для краевых задач 2D точки, 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) непрерывно на интервале интегрирования, заданного в начальной mesh solinit.x и имеет непрерывную первую производную там. Можно использовать sol с deval функция, чтобы оценить решение в других точках в интервале.

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

Поле Описание

x

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

y

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

yp

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

parameters

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

solver

'bvp4c'

stats

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

Больше о

свернуть все

Многоточечные Краевые задачи

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

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

yp = odefun(x,y,k) 

В граничных условиях функционируют bcfun(yleft,yright), yleft(:,k) решение на левом контуре [a k −1, a k]. Точно так же 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 код конечной разности, который реализует трехэтапную формулу [1], [2] Lobatto IIIa. Это - формула словосочетания, и полином словосочетания предоставляет решение C1-continuous, которое является четвертым порядком, точным однородно в интервале интегрирования. Поймайте в сети выбор, и контроль ошибок основаны на невязке непрерывного решения.

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

Ссылки

[1] Шемпин, L.F., и Дж. Кирженка. "Решатель BVP на основе остаточного управления и PSE MATLAB". Математика Сделки ACM. Softw. Издание 27, Номер 3, 2001, стр 299–316.

[2] Шемпин, L.F., М.В. Рейчелт и Дж. Кирженка. "Решая Краевые задачи для Обыкновенных дифференциальных уравнений в MATLAB с bvp4c". MATLAB File Exchange, 2004.

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