Вывод ступенчатой ​​функции Хевисайда в качестве параметра для другой функции Matlab

Я пытаюсь реализовать математическую модель, описанную в этой ссылке:

Чад Лю, Чуан-Юань Ли, Фан Юань. Математическое моделирование восходящего пути Феникса. Plos One, Computational Biology.v.10, 2014.

Я хотел бы, чтобы выходная переменная dPLAdt вычислялась численно, подчиняясь функции Хевисайда, которая вычисляет переменные C3 и C7. Как видно из приведенного ниже кода, функция Хевисайда будет активирована по истечении времени 1440. Программа работает, но не выдает результатов, соответствующих прикрепленной картинке из статьи![http://i.stack.imgur.com/37Vjo.jpg][1]. Я хотел бы сопоставить результат бумаги.

Вот мой код:

1) Функция

function dPLAdt = PLAprod(t,PLA)


dPLAdt = zeros(2,1);

k1 = 2.8*10^(-5);
k3 = 3.24*10^(-6);

kminus5 = 0.09/60;
k5 = 0.06/60;
k2short = 144/60;
k2 = 11;
k4short = 26/60;
k4 = 12;
k7 = 0.06;
p = [k1*heaviside(t-1440); k3*heaviside(t-1440)]; % CASP 3 e CASP7
C3 = p(1); 
C7 = p(2);



dPLAdt(1) = kminus5 - k5*(PLA(1)) - (k2short*(PLA(1))*(C3))/(k2 + PLA(1)) - (k4short*(PLA(1))*(C7))/(k4 + PLA(1));


dPLAdt(2) = (k2short*(PLA(1))*(C3))/(k2 + (PLA(1))) + (k4short*(PLA(1))*(C7))/(k4 + PLA(1)) - k7*(PLA(2));

end

2) Скрипт для вызова решателя

[T,Y] = ode45(@PLAprod,[0,2880],[1500 500]);
plot(T,Y(:,1),'-r', T,Y(:,2),'-b')
legend('iPLA','aPLA')

matlab ode numerical-integration
person Alexandre Sarmento    schedule 26.08.2014    source источник

Ответы (1)


arrow_upward
1
arrow_downward

функция Хевисайда эффективно делает это (очень) жесткая система, которая ode45 не предназначен для обработки. Вы можете попытаться настроить допуски интеграции или использовать ode15s или другой жесткий решатель. , но на самом деле вы имеете дело с кусочно-временной системой. Вы можете точно и эффективно интегрировать его с помощью общего адаптивного решателя, такого как ode45, если разбить время (tspan) и создать две немного разные функции интеграции. Просто замените функцию Хевисайда постоянным низким или высоким значением. Затем смоделируйте первый от 0 до 1440, а второй от 1440 до 2880, используя в качестве начальных условий конечное состояние первого. См. мой ответ на этот вопрос для примера.

Вы можете думать о функции Хевисайда как об операторе if, и почти всегда плохая идея помещать if или любую другую форму ветвления внутрь функции интеграции. Решатели ОДУ предполагают определенную степень гладкости интегрируемых функций.

person horchler    schedule 26.08.2014