Изменить константу в вычислениях ОДУ при определенных условиях с помощью флага

У меня есть ODE для расчета изменения кислотности. Все работает нормально, только хотелось бы менять константу всякий раз, когда кислотность достигает критической точки. Предполагается, что это какой-то необратимый эффект, который я хочу смоделировать.

Мои константы берутся из файла структуры (c), который я загружаю один раз в функции ODE.

[Time,Results] = ode15s(@(x, c) f1(x, c),[0 c.length],x0,options);

Основная проблема, с которой я столкнулся, заключается не в том, чтобы сказать Matlab изменить константу, а в том, чтобы запомнить, если это уже произошло во время симуляции один раз. поэтому Matlab должен взять необратимо измененную константу, а не ту, которую я указал в начале.

Я хотел бы написать флаг, который сохраняется во время работы ODE, и условие if «если флаг существует, измените константу». Как это сделать?

ОБНОВЛЕНИЕ I: ПРОБЛЕМА РЕШЕНА

Здесь первое собственное решение, оно не отшлифовано и требует подхода к структурному файлу. Это означает, что константы, которые должны внезапно измениться при событии, должны быть структурными файлами, которые передаются функцией ODE в функцию, которая должна быть оценена (см. выше синтаксис ODE). Функция принимает такие входные данные:

function [OUTPUT] = f1(t, x, c)

% So here, the constants all start with c. followed by the variable name. (because they are structs instead of globals)

%% write a flag when that event happens

if c.someODEevent <= 999 && exist ('flag.txt')  == 0
    dlmwrite ('flag.txt',1);
end

%% next iteration will either write the flag again or not. more importantly, if it exists, the constant will change because of this.

if exist ('flag.txt')  == 2
        c.changingconstant = c.changingconstant/2;
end



end

Пожалуйста, ознакомьтесь с любезным ответом Horchler, где вы должны позаботиться о том, чтобы такой шаг мог привести к неточностям, и вы должны быть осторожны, чтобы проверить, делает ли ваш код то, что он должен делать.


matlab simulation ode numerical-integration
person Easyquestionsonly    schedule 10.02.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Чтобы сделать это точно, вы должны использовать обнаружение событий в решателе ODE. Я не могу дать вам конкретный ответ, потому что вы указали только ode15s вызов в своем вопросе, но вам нужно будет написать функцию events, а затем указать ее через odeset. Что-то вроде этого:

function acidity_main
% load c data
...
x0 = ...
options = odeset('Events',@events); % add any other options too

% integrate until critical value and stop
[Time1,Results1] = ode15s(@(x,c)f1(x,c),[0 c.length],x0,options);

x0 = Results(end,:); % set new initial conditions
% pass new parameters -it's not clear how you're passing parameters
% if desired, change options to turn off events for faster integration
[Time2,Results2] = ode15s(@(x,c)f1(x,c),[0 c.length],x0,options);

% append outputs, last of 1 is same as first of 2
Time = [Time1;Time2(2:end)];
Results = [Results1;Results2(2:end,:)];
...

function y=f1(x,c)
% your integration function
...

function [value,isterminal,direction] = events(x,c)
value = ... % crossing condition, evaluates to zero at event condition
isterminal = 1; % stop integration when event detected
direction = ... % see documentation

Вы захотите использовать события для интеграции прямо до точки, где «кислотность достигает критической точки», и остановить интеграцию. Затем снова вызовите ode15s с новым значением и продолжите интеграцию. Это может показаться грубым, но это то, как такие вещи можно сделать точно.

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

Примечание: мне кажется странным, что вы передаете то, что вы называете "константами", c, в качестве второго аргумента в ode15s. Эта функция имеет строгие требования к входным аргументам: первая — независимая переменная (часто время), а вторая — массив переменных состояния (такой же, как ваш начальный вектор условия). Также, если f1 принимает только два аргумента, @(x,c)f1(x,c) лишний — достаточно передать @f1.

person horchler    schedule 10.02.2015
comment
Спасибо, я боялся, что мне придется вызывать две функции ODE. Я пытался избежать этого, хотя, нет ли более простого способа? Я стараюсь не усложнять и без того сложный сценарий. - person Easyquestionsonly; 11.02.2015
comment
Передача констант необходима, потому что мне нужно было сделать константы доступными для редактирования на листе Excel, импортировать как переменные в рабочее пространство Matlab и, наконец, превратить их в структуры для передачи в функцию. Вся причина заключалась в том, чтобы избежать использования Matlab при редактировании скрипта и избежать использования глобальных переменных (эти константы должны были быть известны как минимум двум функциям), которые потребляли много памяти. См. здесь: stackoverflow.com/questions/27218341/ - person Easyquestionsonly; 11.02.2015
comment
@Easyquestionsonly: другого пути нет, если только вы не хотите получить неточные и, возможно, совершенно неправильные результаты. Да, вы не должны использовать глобальные. Ссылка, которую вы дали, предназначена для fsolve — целевая функция имеет только один обязательный вход. Решатели ODE, такие как ode15s, имеют два. Если c является фиксированным параметром или их набором, то они должны быть переданы в качестве третьего аргумента. См. этот мой ответ, где n передается в качестве параметра функции интеграции f. - person horchler; 11.02.2015
comment
Мне удалось найти другой способ сделать это. Он имеет общий вид, в котором функция, которую загружает ODE, имеет в начале команду load c.mat. это с. struct содержит мои переменные. Как только я передаю критическое значение в функцию, я сохраняю c.mat, иначе я ничего не делаю. Когда функция снова вызывается в том же цикле, первое, что она видит, это загрузка c.mat с измененным значением. Таким образом, у меня есть одно событие, сохраненное в том же цикле ODE. - person Easyquestionsonly; 12.02.2015
comment
@Easyquestionsonly: я этого не понимаю, потому что вы не предоставили никакого кода, но если он работает, и вы можете подтвердить, что ваши результаты точны и правильны: отлично. Ключевым моментом является то, что вы не вводите разрывы, шаги, операторы if и т. д. в свою функцию интеграции, функцию, которая описывает ваши ODE и вызывается ode15s (f1 в случае вашего вопроса). Это сделает вашу систему более жесткой, снизит точность и, возможно, приведет к совершенно неправильным результатам. - person horchler; 12.02.2015
comment
Спасибо, что указали на это. Суть в том, чтобы ввести какое-то возмущение (фактически отреагировать на возмущение) в довольно динамичной среде. Артефакты можно легко распознать в моем случае, но, возможно, не в вашем. - person Easyquestionsonly; 12.02.2015