在Matlab中求解微分方程

Question

我需要同时求解这两个微分方程。

dr^3/dt=(-3*D*Cs)/(ρ*r0^2 )*r*(1-C)

dC/dt=((D*4π*r0*N*(1-C)*r)-(Af*C))/V

注意：dr ^ 3 / dt是r ^ 3对t的导数

这两个方程类似于微粒悬浮液的溶解过程及其在血液中的同时吸收随时间变化的粒子半径（r）和浓度（C）。当固体溶解时，随着半径的减小，半径r将减小，浓度C将增加，最终达到平稳状态（即达到平衡），因为通过此Af * C将溶解的固体移至血液中项（其中Af是某种吸收速率常数）。这些方程式来自本文，我正在尝试复制它们：jpharmsci.org/article/S0022-3549(18)30334-4/fulltext#sec3.2.1-C随t的变化应该类似于图3（DCU示例）。

我做了简化：dr ^ 3 / dt = 3r ^ 2 *（dr / dt），然后将方程式的两边除以3r ^ 2。颂歌变成：

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

y(1) = r* and 
y(2) = C*

r* and C* 

是本文中使用的术语，是“归一化”的半径和浓度，其中

r*=r/r0 and C*=C/Cs

其中：

• r =粒子半径（随时间变化并用dydt（1）表示）
• r0 =初始粒子半径
• C =溶解固体的浓度（随时间变化并用dydt（2）表示）
• Cs =饱和溶解度

其余代码如下。 根据作者对纸上使用的值的反馈进行更新，并将初始值校正为y0 = [1 0]

MW=224; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60] equation provided by authors, divide by 600,000 to convert to m2/s 
rho=1300; %kg/m3
r0=10.1e-6; %m dv50
Cs=1.6*1e6/1e9; %kg/m3 - 1.6ug/m3 converted to kg/m3
V=5*0.3/1e6;%m3 5ml/Kg animal * 0.3Kg animal, divide by 1e6 to convert to m3
W=30*0.3/1000000; %kg; 30mg/Kg animal * 0.3Kg animal, divide by 1e6 to convert to m3
N=W/((4/3)*pi*r0^3*rho); % particle number
Af=0.7e-6/60; %m3/s
tspan=[0 24*3600]; %s in 24 hrs
y0=[1 0];
[t,y]=ode113(@(t,y) odefcnNY_v11(t,y,D,Cs,rho,r0,Af,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*');
plot(t/3600, y(:,2)*Cs) % time in hr, and bulk concentration on y
xlabel('time, hr')
ylabel('C, kg/m3')
legend('Dissolved drug concentration')
title ('C');

我首先尝试使用ode45，但是代码花费了很长时间才能运行，最终出现一些错误。然后，我尝试使用ode113并收到以下错误。

Warning: Failure at t=2.112013e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

更新：已更新功能代码以解决奇异性问题：

function dydt=odefcnNY_v10(t,y,D,Cs,rho,r0,N,V,Af)
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)-Af*y(2))/V; %dC*/dt
end

结果

模型的机械背景

dr / dt的导数

dC / dt的导数

模型假设

以上，您将找到显示这些方程式推导的幻灯片。他们假设在Noyes-Whitney方程中，溶出速率dM / dt =（φ/ h）4 ^^ 2φ（φ-φ），膜厚h等于颗粒半径r。如果雷诺数低且颗粒小于60um（在这种情况下为10um），这通常是在生物制药模型中进行的简化。如果做这个假设，我们剩下dM / dt =𝐷4 𝜋 𝑟𝑁（𝑁−𝐶𝑠）。我渴望复制本文，因为我想做同样的事情，即对皮下注射微悬液的药物吸收进行建模。我已经联系了作者，他们似乎不太确定自己做了什么，因此我正在寻找其他资料来源，例如本文：https://pubs.acs.org/doi/pdf/10.1021/acs.iecr.7b04730在等式6中，显示了dC / dt的等式。他们将单位体积的表面积（a）（等式5）的变化嵌入等式6中。它们的传质系数kL是集总参数= D / h（扩散率/膜厚度）。

我需要同时求解这两个微分方程。 dr ^ 3 / dt =（-3 * D * Cs）/（ρ* r0 ^ 2）* r *（1-C）dC / dt =（（D *4π* r0 * N *（1-C）* r）-（Af * C））/ V注意：dr ^ 3 / dt是r ^ 3相对于t的导数...

Answer 1

以半径方程的原始形式