Я пытаюсь решить проблему, аналогичную этой: Решение дифференциальных уравнений в 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 часа и все еще не завершен. Что теперь отличается от предыдущего вопроса по ссылке выше, что делает это так долго?
Спасибо
odeset
для явного определения абсолютного и относительного допуска. Если я установлю их высокими, решатель действительно займет много времени. Допуски по умолчанию должны быть достаточно низкими, но проверьте это или установите их. Abstol до 1e-6 и RelTol до 1e-9 или меньше. - person Lutz Lehmann   schedule 29.11.2019options = 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