Matlab: возможно ли численно решить систему од со смесью начальных и конечных условий?

Я пытаюсь использовать ode45 для решения системы ОДУ:

[X,Y]=  ode45(@sys,[0, T],y0);

куда,

function dy = sys(t,y)

        dy(1) = f_1(y)
        dy(2) = f_2(y)
        dy(3) = f_3(y)
end

Проблема в том, что функция ode45 требует, чтобы y0 были начальными значениями [y_1(0), y_2(0), y_3(0)], тогда как в моей системе доступны только значения [y_2(0), y_3(0), y_3(T)].

Математически этого набора начальных/конечных условий должно быть достаточно, чтобы зафиксировать систему, но есть ли способ, которым я могу работать с этим с помощью ode45 или любых других функций в MATLAB?

Спасибо!


matlab numerical-methods ode
person Vokram    schedule 14.06.2013    source источник
comment
Меня очень интересует этот вопрос, но, боюсь, я не смогу помочь; Я никогда раньше не сталкивался с такой проблемой... Я знаю, что ode45 можно интегрировать в обратном направлении (просто используйте tspan = [tend tstart]), поэтому вы можете придумать итеративную схему для получения y_1, чтобы y_3(0) и y_3(T) были удовлетворены. Излишне говорить, что это, вероятно, будет очень медленным и довольно неуклюжим, но это будет решение. Я буду следить за этим :) Не могли бы вы опубликовать уравнения и начальные/конечные условия?   -  person Rody Oldenhuis    schedule 14.06.2013
comment
@RodyOldenhuis Думаю, я нашел способ решить эту проблему.   -  person Vokram    schedule 14.06.2013

Ответы (2)


arrow_upward
6
arrow_downward

Немного покопавшись в документации Matlab, я думаю, что более элегантным способом является использование функции bvp4c. bvp4c — это функция, специально предназначенная для решения подобных проблем с граничными значениями, в отличие от ode**, которые на самом деле предназначены только для проблем с начальными значениями. На самом деле в Matlab есть целый набор других функций, таких как deval и bvpinit, которые действительно облегчают использование bvp4c. Вот ссылка на документацию Matlab.

Я опубликую краткий (и, возможно, немного надуманный) пример здесь:

function [y1, y2, y3] = test(start,T)

solinit = bvpinit(linspace(0,3,10), [1,1,0]);
sol = bvp4c(@odefun,@bvpbc,solinit);

tspan = linspace(start,T,100);
S = deval(sol, tspan);
y1 = S(1,:);
y2 = S(2,:);
y3 = S(3,:);

plot (tspan,y1)

figure
plot (tspan,y2)

figure
plot (tspan,y3)


%% system definition & BVCs

function dydx = odefun(t,y)
    dydx(1) = y(1) + y(2) + t;

    dydx(2) = 2*y(1) + y(2);

    dydx(3) = 3 * y(1) - y(2);
end

function res = bvpbc(y0,yT)
   res= [y0(3) yT(2) yT(3)];
end

end

Функция test выводит 3 вектора точек решений для y1, y2 и y3 соответственно, а также строит их.

Вот пути переменных, построенные Matlab: y1 pathпуть y2 путь y3

Кроме того, я нашел очень полезным это видео профессора Джейка Бланшара из WMU. .

person Vokram    schedule 14.06.2013
comment
+2: похоже, ты лизнул его! Сегодня я узнал кое-что новое, спасибо за это :) - person Rody Oldenhuis; 14.06.2013
comment
+2 за удивление, что @RodyOldenhuis не знал о решателе BVP! :-) Вы можете принять свой собственный ответ, если хотите. - person horchler; 15.06.2013

arrow_upward
3
arrow_downward

Один из подходов заключается в использовании метода стрельбы для итеративного поиска неизвестного начального состояния y_1(0) таким образом, чтобы желаемое конечное состояние y_3(T) удовлетворено.

Итерация продолжается путем решения нелинейного уравнения F = 0:

F(y_1(0)) = Y_3(T) - y_3(T)

где функция Y в верхнем регистре — это решение, полученное путем интегрирования вперед ОДУ по времени из набора начальных условий. Задача состоит в том, чтобы угадать значение y_1(0), чтобы получить F = 0.

Поскольку сейчас мы решаем нелинейное уравнение, применимы все обычные подходы. В частности, вы можете использовать метод деления пополам или метод Ньютона, чтобы обновить свое предположение для неизвестного начального состояния y_1(0). Обратите внимание на пару вещей:

  1. ОДУ интегрируются в [0,T] на каждой итерации нелинейного решения.
  2. Для F = 0 может быть несколько решений, в зависимости от структуры вашего ODE.
  3. Метод Ньютона может сходиться быстрее, чем пополам, но также может быть численно нестабильным, если вы не можете дать хорошее начальное предположение для y_1(0).

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

Надеюсь это поможет.

person Darren Engwirda    schedule 14.06.2013
comment
+1: это был бы мой подход, но я думаю, что собственный ответ @Vokram должен быть проголосован за первое место :) - person Rody Oldenhuis; 14.06.2013