Этот пример показывает, как решить численно трудную краевую задачу с помощью продолжения, которое эффективно разбивает проблему в последовательность более простых проблем.
Для , рассмотрите дифференциальное уравнение
.
Проблема создана на интервале и подчиняется граничным условиям
,
.
Когда , решение уравнения подвергается быстрому переходу рядом , таким образом, трудно решить численно. Вместо этого этот пример использует продолжение, чтобы выполнить итерации через несколько значений до . Промежуточные решения каждый используются в качестве исходного предположения для следующей проблемы.
Чтобы решить эту систему уравнений в MATLAB, необходимо закодировать уравнения, граничные условия и исходное предположение прежде, чем вызвать решатель для краевой задачи bvp4c
. Вы любой может включать необходимые функции как локальные функции в конце файла (как сделано здесь), или можно сохранить их как отдельные, именованные файлы в директории на пути MATLAB.
Используя замены и , можно переписать уравнение как систему уравнений первого порядка
,
.
Запишите функцию, чтобы закодировать уравнения с подписью dydx = shockode(x,y)
, где:
x
является независимой переменной.
y
является зависимой переменной.
dydx(1)
дает уравнение для , и dydx(2)
дает уравнение для .
Сделайте функцию векторизованной, так, чтобы shockode([x1 x2 ...],[y1 y2 ...])
возвратил [shockode(x1,y1) shockode(x2,y2) ...]
. Этот подход улучшает производительность решателя.
Соответствующая функция
function dydx = shockode(x,y) pix = pi*x; dydx = [y(2,:) -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)]; end
Примечание: Все функции включены как локальные функции в конце примера.
Решатель BVP требует, чтобы граничные условия были в форме . В этой форме граничные условия:
,
.
Запишите функцию, чтобы закодировать граничные условия с подписью res = shockbc(ya,yb)
, где:
ya
является значением граничного условия в начале интервала .
yb
является значением граничного условия в конце интервала .
Соответствующая функция
function res = shockbc(ya,yb) % boundary conditions res = [ya(1)+2 yb(1)]; end
Аналитические Якобианы для функции ОДУ и граничных условий могут быть вычислены легко в этой проблеме. Обеспечение Якобианов делает решатель более эффективным, поскольку решатель больше не должен аппроксимировать их с конечными разностями.
Для функции ОДУ якобиевская матрица
.
Соответствующая функция
function jac = shockjac(x,y,e) jac = [0 1 0 -x/e]; end
Точно так же для граничных условий, якобиевские матрицы
, .
Соответствующая функция
function [dBCdya,dBCdyb] = shockbcjac(ya,yb) dBCdya = [1 0; 0 0]; dBCdyb = [0 0; 1 0]; end
Используйте постоянное предположение для решения на сетке пяти точек в .
sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);
При попытке решить уравнение непосредственно с помощью , затем решатель не может преодолеть плохое создание условий проблемы около точка перехода. Вместо этого чтобы получить решение для , этот пример использует продолжение путем решения последовательности проблем для , , и . Вывод от решателя в каждой итерации действия как предположение для решения в следующей итерации (это - то, почему переменной для исходного предположения от bvpinit
является sol
и вывод от решателя, также называют sol
).
Поскольку значение якобиана зависит от значения , установите опции в цикле, задающем shockjac
и функции shockbcjac
для Якобианов. Кроме того, включите векторизацию, поскольку shockode
закодирован, чтобы обработать векторы значений.
e = 0.1; for i = 2:4 e = e/10; options = bvpset('FJacobian',@(x,y) shockjac(x,y,e),'BCJacobian',@shockbcjac,'Vectorized','on'); sol = bvp4c(@(x,y) shockode(x,y,e),@shockbc, sol, options); end
Постройте результаты bvp4c
для mesh и решение . С продолжением решатель может обработать разрыв в .
plot(sol.x,sol.y(1,:),'-o'); axis([-1 1 -2.2 2.2]); title(['There Is a Shock at x = 0 When e =' sprintf('%.e',e) '.']); xlabel('x'); ylabel('solution y');
Перечисленный здесь локальные функции, которые решатель BVP bvp4c
вызывает, чтобы вычислить решение. Также можно сохранить эти функции как их собственные файлы в директории на пути MATLAB.
function dydx = shockode(x,y,e) % equation to solve pix = pi*x; dydx = [y(2,:) -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)]; end %------------------------------------------- function res = shockbc(ya,yb) % boundary conditions res = [ya(1)+2 yb(1)]; end %------------------------------------------- function jac = shockjac(x,y,e) % jacobian of shockode jac = [0 1 0 -x/e]; end %------------------------------------------- function [dBCdya,dBCdyb] = shockbcjac(ya,yb) % jacobian of shockbc dBCdya = [1 0; 0 0]; dBCdyb = [0 0; 1 0]; end %-------------------------------------------