Решение дифференциальных уравнений в Matlab — In-Vitro Dissolution

Я пытаюсь решить проблему, аналогичную этой: Решение дифференциальных уравнений в Matlab

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

function dydt=odefcnNY_v12(t,y,D,Cs,rho,r0,N,V)
dydt=zeros(2,1);
dydt(1)=(-D*Cs)/(rho*r0^2)*(1-y(2))*y(1)/(1e-6+y(1)^2); % dr*/dt
dydt(2)=(D*4*pi*N*r0*(1-y(2))*y(1))/V; %dC*/dt
end

т. е. член поглощения из предыдущего вопроса удаляется:

Absorption term: Af*y(2)

Соединение также отличается, поэтому MW, Cs и r0 отличаются, и экспериментальная установка также отличается, поэтому W и V теперь изменены. Чтобы учесть эти изменения, вызов ode113 меняется на это:

MW=336.43; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s 
rho=1300; %kg/m3
r0=9.75e-8; %m dv50
Cs=0.032; %kg/m3 
V=0.0009;%m3 900 mL dissolution bath
W=18e-6; %kg 18mg
N=W/((4/3)*pi*r0^3*rho); % particle number
tspan=[0 200*3600]; %s in 200 hours
y0=[1 0];
[t,y]=ode113(@(t,y) odefcnNY_v12(t,y,D,Cs,rho,r0,N,V), tspan, y0);
plot(t/3600,y(:,1),'-o') %plot time in hr, and r*
xlabel('time, hr')
ylabel('r*, (rp/r0)')
legend('DCU')
title ('r*');
plot(t/3600,y(:,1)*r0*1e6); %plot r in microns
xlabel('time, hr');
ylabel('r, microns');
legend('DCU');
title('r');
plot(t/3600,y(:,2),'-') %plot time in hr, and C*
xlabel('time, hr')
ylabel('C* (C/Cs)')
legend('DCU')
title('C*');

Текущая проблема в том, что этот код работает уже 3 часа и все еще не завершен. Что теперь отличается от предыдущего вопроса по ссылке выше, что делает это так долго?

Спасибо

введите здесь описание изображения


person Engineer_1331    schedule 29.11.2019    source источник
comment
Я использовал предыдущий код, установив Af равным нулю в определении новых констант. Нет проблем в интеграции, используя неявный метод Радау. Полное растворение и стабилизация при С*=0,625 в течение первой секунды. Вы действительно должны использовать odeset для явного определения абсолютного и относительного допуска. Если я установлю их высокими, решатель действительно займет много времени. Допуски по умолчанию должны быть достаточно низкими, но проверьте это или установите их. Abstol до 1e-6 и RelTol до 1e-9 или меньше.   -  person Lutz Lehmann    schedule 29.11.2019
comment
Я изменил свой код на этот options = odeset('RelTol',1e-6,'AbsTol',1e-9); [t,y]=ode113(@(t,y) odefcnNY_v11(t,y,D,Cs,rho,r0,N,V,Af), tspan, y0, options);, но его выполнение все еще занимает много времени. Я поставил его на паузу, чтобы посмотреть на ошибки, я прилагаю скриншот в вопросе.   -  person Engineer_1331    schedule 29.11.2019
comment
Пробовали ли вы что-либо из: а) уменьшить интервал интегрирования до 1 секунды, б) перейти на другой решатель, в) ввести глобальную переменную счетчика и распечатать время для каждого сотого вызова функции ОДУ и 10 последующих (или выполнить полный журнал), чтобы получить отчет о ходе работы и выяснить, где интегратор останавливается, или другие меры по отладке?   -  person Lutz Lehmann    schedule 29.11.2019
comment
Мы, наверное, ужасно не по теме. Ваш код работает без ошибок, просто результаты не такие, как ожидалось. Или что численный решатель останавливается для этого конкретного ODE. Что больше связано с ODE, чем с решателем. Такие обсуждения лучше размещать на scicomp.SE.   -  person Lutz Lehmann    schedule 29.11.2019


Ответы (1)


Я не могу воспроизвести вашу проблему. Я использую "стандартные" модули python numpy и scipy, скопировал блок параметров,

MW=336.43; # molecular weight
D=9.916e-5*(MW**-0.4569)*60/600000 #m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s 
rho=1300.; #kg/m3
r0=9.75e-8; #m dv50
Cs=0.032; #kg/m3 
V=0.0009;#m3 900 mL dissolution bath
W=18e-6; #kg 18mg
N=W/((4./3)*pi*r0**3*rho); # particle number
Af = 0; # bath is isolated 

использовал ту же функцию ODE, что и в предыдущем посте (помните Af=0)

def odefcnNY(t,y,D,Cs,rho,r0,N,V,Af):
    r,C = y;
    drdt = (-D*Cs)/(rho*r0**2)*(1-C) * r/(1e-10+r**2); # dr*/dt
    dCdt = (D*4*pi*N*r0*(1-C)*r-(Af*C))/V;            # dC*/dt
    return [ drdt, dCdt ];

и решил ОДУ

tspan=[0, 1.0]; #1 sec
#tspan=[0, 200*3600]; #s in 200 hours
y0=[1.0, 0.0];
method="Radau"

sol=solve_ivp(lambda t,y: odefcnNY(t,y,D,Cs,rho,r0,N,V,Af), tspan, y0, method=method, atol=1e-8, rtol=1e-11);

t = sol.t; r,C = sol.y;
print(sol.message)
print("C*=",C[-1])

Это работает мгновенно, используя 235 шагов в первую секунду и 6 дополнительных шагов, чтобы охватить постоянное поведение в оставшееся время из 200 часов.

графики компонентов

Я могу только испортить это, увеличив допуски до необоснованно больших значений, таких как 1e-4, и только если эпсилон, используемый в смягчении, равен 1e-12. Тогда жесткий поворот, когда радиус достигает нуля, становится слишком жестким, регулятор размера шага зацикливается. Это скорее ошибка грубой реализации контроллера размера шага, чего не должно быть в подпрограммах Matlab.

person Lutz Lehmann    schedule 29.11.2019
comment
Ты прав. Я переключился на ode15s, и это решило проблему. - person Engineer_1331; 29.11.2019