Решите BVP Используя продолжение

Этот пример показывает, как решить численно трудную краевую задачу с помощью продолжения, которое эффективно разбивает проблему в последовательность более простых проблем.

Для 0<e1, рассмотрите дифференциальное уравнение

ey+x, y =-eπ2потому что(πx)-πxsin(πx).

Проблема создана на интервале [-1,1] и подчиняется граничным условиям

y(-1)=-2,

y(1)=0.

Когда e=10-4, решение уравнения подвергается быстрому переходу рядом x=0, таким образом, трудно решить численно. Вместо этого этот пример использует продолжение, чтобы выполнить итерации через несколько значений e до e=10-4. Промежуточные решения каждый используются в качестве исходного предположения для следующей проблемы.

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

Уравнение кода

Используя замены y1=y и y2=y, можно переписать уравнение как систему уравнений первого порядка

y1=y2,

y2=-xey-π2потому что(πx)-πxesin(πx).

Запишите функцию, чтобы закодировать уравнения с подписью dydx = shockode(x,y), где:

  • x является независимой переменной.

  • y является зависимой переменной.

  • dydx(1) дает уравнение для y1, и dydx(2) дает уравнение для y2.

Сделайте функцию векторизованной, так, чтобы 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 требует, чтобы граничные условия были в форме g(y(a),y(b))=0. В этой форме граничные условия:

y(-1)+2=0,

y(1)=0.

Запишите функцию, чтобы закодировать граничные условия с подписью res = shockbc(ya,yb), где:

  • ya является значением граничного условия в начале интервала [a,b].

  • yb является значением граничного условия в конце интервала [a,b].

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

function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end

Якобианы кода

Аналитические Якобианы для функции ОДУ и граничных условий могут быть вычислены легко в этой проблеме. Обеспечение Якобианов делает решатель более эффективным, поскольку решатель больше не должен аппроксимировать их с конечными разностями.

Для функции ОДУ якобиевская матрица

JОДУ=fy=[f1y1f1y2f2y1f2y2]=[010-xe].

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

function jac = shockjac(x,y,e)
jac = [0   1
       0  -x/e];
end

Точно так же для граничных условий, якобиевские матрицы

Jy(a)=[1000] , Jy(b)=[0010].

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

function [dBCdya,dBCdyb] = shockbcjac(ya,yb)
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end

Сформируйте исходное предположение

Используйте постоянное предположение для решения на сетке пяти точек в [-1,1].

sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);

Решите уравнение

При попытке решить уравнение непосредственно с помощью e=10-4, затем решатель не может преодолеть плохое создание условий проблемы около x=0 точка перехода. Вместо этого чтобы получить решение для e=10-4, этот пример использует продолжение путем решения последовательности проблем для 10-2, 10-3, и 10-4. Вывод от решателя в каждой итерации действия как предположение для решения в следующей итерации (это - то, почему переменной для исходного предположения от bvpinit является sol и вывод от решателя, также называют sol).

Поскольку значение якобиана зависит от значения e, установите опции в цикле, задающем 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 x и решение y(x). С продолжением решатель может обработать разрыв в x=0.

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');

Localfunctions

Перечисленный здесь локальные функции, которые решатель 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
%-------------------------------------------

Смотрите также

| |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте