在 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*

是论文中使用的术语，是 "normalised" 半径和浓度，其中

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（diffusivity/film厚度）。

Answer 1

半径方程的原始形式

d(r^3)/dt = -3K*(r^3)^(1/3)*(1-C)

或低功率的

dr/dt = -K/r*(1-C)  <==> d(r^2)/dt = -2K*(1-C)

你在半径缩小到零的那一刻达到了一个奇点，就像大约r=sqrt(A-B*t)，也就是说，小球会在有限的时间内消失。从那一刻起，显然，一个人应该有 r=0。这可以通过将模型从 2 组件系统更改为仅 C 的标量方程来获得，通过事件机制获得准确的时间。

或者，您可以修改半径方程，使 r=0 成为自然静止点。第一个版本以某种方式做到了这一点，但由于第二个版本是等效的，因此仍然可以预料到数字上的困难。可能的修改是

d(r^2)/dt = -2K*sign(r)*(1-C)

其中 signum 函数可以用连续近似代替，例如

x/(eps+abs(x)), x/max(eps,abs(x)), tanh(x/eps), ...

或者在简化形式中，奇点可以缓和为

dr/dt = -K*(1-C) * r/(eps^2+r^2)

或仍在其他变体中

dr/dt = -K*(1-C) * 2*r/(max(eps^2,r^2)+r^2)

应用于给出的具体案例（使用python代替matlab）

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

并鉴于 r=0 处的近奇点应用隐式方法

D=8.3658e-10#m2/s 
rho=1300; #kg/m3
r0=10.1e-6; #m dv50
Cs=0.0016; #kg/m3
V=1.5e-6;#m3
W=9e-6; #kg
N=W/(4/3*pi*r0^3*rho);
Af=0.7e-6/60; #m3/s
tspan=[0, 24*3600]; #sec in 24 hours
y0=[1.0, 0.0]; # relative radius starts at full, 1.0
sol=solve_ivp(lambda t,y: odefcnNY(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0, method="Radau", atol=1e-14);
t = sol.t; r,C = sol.y;

然后生成解图

经过修正的参数现在看起来与发布的图表很接近。

在 Matlab 中求解微分方程

Solving differential equations in Matlab

matlab

ode

ode45