ode45 сходится к правильной форме кривой, но с неправильным решением - PullRequest
Новая
       104

ode45 сходится к правильной форме кривой, но с неправильным решением

0 голосов
/ 15 февраля 2019

Заранее спасибо за помощь.Я не ищу явного решения своей проблемы, а скорее указал на мои, вероятно, очевидные ошибки.

Я подключился к решению системы нелинейных ОДУ первого порядка в MATLAB.Система была решена численно в этом исследовании: http://web.math.ku.dk/~moller/e04/bio/ludwig78.pdf

Я следовал документации для ode45 , и у меня есть код, который выполняется.

Я сделал всеработа, чтобы понять и воссоздать модель с нуля.Я представил качественную часть для классного проекта.То, что я делаю сейчас, - это продвижение этого проекта на шаг дальше, решая систему в MATLAB с помощью runge-kutta (или любого другого метода, который работает).Наконец, я хочу углубиться в теорию численного анализа, чтобы выяснить, почему выбранный метод сходится.

Вот график численно решенной системы, которую я пытаюсь воссоздать: enter image description here

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

  1. Шкала времени, в течение которой происходит изменение, в три раза больше, чем на приведенном выше графике.
  2. Диапазон значений функцииэто совершенно неправильно.
  3. Желаемые формы появляются, только если я настраиваю начальные условия, чтобы они значительно отличались от того, что показано около t = 0 выше.

Так что я смотрюявляется причиной этих расхождений.Я столько раз проверял свою систему ODE и значений параметров, что мои глаза размыты.Возможно, я что-то упускаю концептуально?

Код: 

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

tspan = [0 200];
init = [1 1 1];

[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);
subplot(3,1,1);
plot(t,Y(:,1),'b');
title('Budworm Density');

subplot(3,1,2)
plot(t,Y(:,2),'g');
title('Branch Density');

subplot(3,1,3);
plot(t,Y(:,3),'r');
title('Foliage Condition');

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)
dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
         r_s*y(2)*(1 - (y(2)*k_e)/(k_s*y(3)));
         r_e*y(3)*(1 - (y(3)/k_e)) - p*y(1)/y(2)
        ];

end

1 Ответ

0 голосов
/ 20 февраля 2019

Я не вижу ничего плохого в вашем коде как таковом.Но я думаю, что при создании рисунка есть некоторые тонкости, которые недостаточно хорошо объяснены в статье.

1) Ось S масштабируется (на этикетке написано «относительный»).Я считаю, что они масштабировали S на k_s.Я думаю, что вам также нужно масштабировать параметр p (установить p = p * k_s), иначе последний член в уравнении для E будет крошечным, и популяция E не уменьшится в течение требуемых временных масштабов.

2Я думаю, что они должны были установить некоторый нижний предел для E, чтобы избежать деления на 0. Вы можете видеть на рисунке, что сначала E-> 0, но в вашем уравнении для S, если бы это произошло, вы бы делили на 0 ирешатель не сходится.

Собрав их вместе, следующая небольшая модификация вашего кода даст результат, более похожий на тот, который показан в статье:

ODE result

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

% Scale p with k_s
p = p*k_s;

tspan = [0 50]; % [0 200];
init = [1e-16 0.075*k_s 1]; % [1 1 1];

[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);

% To scale before plotting, so everything fits on a 0->1 y axis.
maxB = 500;
S_scale = k_s;

figure('Position', [200 200 1000 600]);
hold on;
plot(t,Y(:,1)/maxB,'b');
plot(t,Y(:,2)/(S_scale),'g');
plot(t,Y(:,3),'r');

ylim([0, 1]);

hold off;

box on;

legend({['Budworm Density, B / ', num2str(maxB)], 'Branch Density, S / 0.75', 'Foliage Condition, E'}, ...
    'Location', 'eastoutside')

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)

% Place lower limit on E
E = max(y(3), 1e-5);

dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
         r_s*y(2)*(1 - (y(2)*k_e)/(k_s*E));
         r_e*E*(1 - (E/k_e)) - p*y(1)/y(2)
        ];

end

Существует большая чувствительность к начальным условиям.

Дальнейшая настройка еще ближе приближает вас к исходной фигуре, но я не уверен, что это всего лишь взлом: в первом уравнении замените k_b * y (2) на просто k_b.Без этого плотность Будды становится слишком большой, прежде чем уменьшится.Новый участок ниже.

Modified k_b

...