Обновить переменные ODE в Matlab

Question

Я пытаюсь реплицировать результаты статьи (Immune Response). Короче говоря, он рассматривает концентрацию патогенов, антител и плазматических клеток, а также влияние иммунного ответа на орган. В зависимости от состояния органа изменяется переменная, которая влияет на концентрацию плазматических клеток.

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

Как обновить значения, чтобы они интерпретировались ode45 (или ode23) и есть ли лучший подход для решения этой проблемы?

Мой код ниже.

function dx = iir_2(t, x) 
dx = [0; 0; 0; 0]; 
a11 = 1; 
a31 = 1; 
a41 = 1; 
a12 = 1; 
a22 = 3; 
a32 = 1.5; 
a42 = 1; 
a23 = 1; 
a33 = 0.5; 

b1 = -1; 
b2 = 1; 
b3 = 1; 
b4 = -1; 

u1 = 0; 
u2 = 0; 
u3 = 0; 
u4 = 0; 

tau = 0.1; 

if x(4) >= 0.5 
    a21x4 = 0; 
else 
    a21x4 = cos(pi * x(4)); 
end 

if x(4) > 1 
    x(4) = 1; 
end 

dx(1) = (a11 - a12 * x(3)) * x(1) + b1 * u1; 
dx(2) = a21x4 * a22 * x(1) * (t - tau) * x(3) * (t - tau) - a23 * (x(2) - 2) + b2 * u2; 
dx(3) = a31 * x(2) - (a32 + a33 * x(1)) * x(3) + b3 * u3; 
dx(4) = a41 * x(1) - a42 * x(4) + b4 * u4; 

Этот второй сценарий, который вызывает выше функция ..

Case = ['Subclinical' 'Clinical' 'Chronic' 'Lethal']; 
x1_0 = [1.5 2 2.57 3]; 
x2_0 = 2; 
x3_0 = (1*x2_0)/1.5; 
x4_0 = 0; 

for i = 1:4 
    if i == 1 
     state = 'Subclinical'; 
    elseif i == 2 
     state = 'Clinical'; 
    elseif i == 3 
     state = 'Chronic'; 
    else 
     state = 'Lethal'; 
    end 

    options = odeset('RelTol', 1e-3, 'NonNegative', [1 2 3 4]); 
    [t,x] = ode45(@iir_2, [0 10], [x1_0(i) x2_0 x3_0 x4_0], options); % use ode23??? 
    figure 
    plot(t,x); 
    str = sprintf('Case Number = %d\nx(0) = %d\n%s', i, x1_0(i), state); 
    title(str); 
    axis([0,10,0,10]) 
    legend('Pathogen', 'Plasma Cell', 'Antibody', 'Organ'); 
end 

условные утверждения в функции.

if x(4) >= 0.5 
     a21x4 = 0; 
    else 
     a21x4 = cos(pi * x(4)); 
    end 

    if x(4) > 1 
     x(4) = 1; 
    end 

Answer 1

Существует два способа сделать это.

В любом случае вам нужно понять, что вы не можете просто установить x (4) в 1 и надеяться, что все будет работать наилучшим образом. Matlab не заботится о предыдущих значениях x (4), поскольку они все хранятся в памяти. Кроме того, предстоящее значение x (4) определяется dx (4) и ранее сохраненными значениями x (4) (которые вы не можете установить).

У вас есть два возможных решения вашей проблемы:

A) множество дх (4) = 0, если ваше условие выполнено

if x(4) > 1 
    dx(4) = 0; 
else 
    dx(4) = a41 * x(1) - a42 * x(4) + b4 * u4; 
end 

Это, однако, не приведет к совершенной х (4) = 1, скорее будет небольшая ошибка.

B) Вызов функции события, которую можно узнать о here и here, и в этом случае вы можете прервать вызов новой функции.