求解Matlab中的微分方程，相对于t具有一个方程y ^ 3

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 *并且y（2）= C *（r *和C *是本文中使用的术语，是“归一化”的半径和浓度，其中r * = rp / r0和C * = C / Cs，其中：

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

其余代码如下。

D=4e-9;%m2/s
rho=1300; %kg/m3
r0=10.1e-6; %m dv50
Cs=0.0016; %kg/m3
V=1.5e-6;%m3
W=4.5e-8; %kg
N=W/(4/3*pi*r0^3*rho);
Af=0.7e-6/60; %m3/s
tspan=[0 24*3600]; %s in 24 hrs
y0=[r0 0];
[t,y]=ode45(@(t,y) odefcnNY_v3(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*')
legend('DCU')
plot(t/3600,y(:,2),'-') %plot time in hr, and C*
xlabel('time,hr')
ylabel('C* (C/Cs)')
legend('DCU')

我首先尝试使用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.

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)

或可能已更正

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

或

dr/dt = -K*(1-C)

您在半径缩小为零的那一刻达到奇点，即小球消失了。从那时起，显然应该有r=0。这可以通过从2分量系统到仅C的标量方程的模型更改来获得，并通过事件机制获得准确的时间。

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

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

或者如果我的纠正是正确的

dr/dt = -K*sign(r)*(1-C)

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

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