在Matlab中求解微分方程-体外溶出度 [英] Solving differential equations in Matlab - In-Vitro Dissolution

问题描述

我正在尝试解决与此问题类似的问题:在Matlab中求解微分方程

I am trying to solve a similar problem to this one: Solving Differential Equations in Matlab

但是，这次的情况不是将药物注射到皮下组织中并随后溶解，而是更简单的情况，允许将悬浮液溶解在900毫升的溶出浴中.

However, this time the scenario is not injection of a drug into the subcutaneous tissue and its subsequent dissolution, but a more simple situation where the suspension is allowed to dissolve in a dissolution bath of volume 900 ml.

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

即上一个问题的吸收项已删除:

i.e. the absorption term from the previous question is removed:

Absorption term: Af*y(2)

化合物也不同，因此MW，Cs和r0也不同，并且实验设置也不同，因此现在更改了W和V.为了允许这些更改，ode113调用对此进行了更改:

The compound is also different, so the MW, Cs and r0 are different, and the experimental setup is also different so W and V are now changed. To allow for these changes, the ode113 call changes to this:

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个小时，但仍未完成.现在与上面链接中的上一个问题有什么不同，就是要花这么长时间?

The current problem is that this code has been running for 3 hours, and still not complete. What is different now to the previous question in the link above, that is making it take so long?

谢谢

推荐答案

我真的无法重现您的问题.我使用标准" python模块numpy和scipy，复制了参数块，

I can not really reproduce your problem. I use the "standard" python modules numpy and scipy, copied the block of parameters,

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)

used the same ODE function like in the previous post (remember 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 ];

并解决了ODE

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个步骤，并在接下来的200个小时中使用了另外6个步骤来覆盖恒定的行为.

This works in a snap, using 235 steps for the first second and 6 further steps to cover the constant behavior in the remaining time of the 200 hours.

我只能通过将容差提高到不合理的大值(例如1e-4)来搞砸，而且前提是用于变现的epsilon为1e-12.然后，当半径达到零时的硬转弯太难了，步长控制器陷入循环.这更多是步长控制器的粗略实现的错误，在Matlab例程中不应该这样.

I can only mess it up by increasing the tolerances to unreasonably large values like 1e-4, and only if the epsilon used in the mollification is 1e-12. Then the hard turn when the radius reaches zero is too hard, the step size controller falls into a loop. This is more an error of the crude implementation of the step size controller, that should not be the case in the Matlab routines.

这篇关于在Matlab中求解微分方程-体外溶出度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！