Matlab求解体外溶出微分方程耗时过长问题咨询

Let's walk through the key differences between your current setup and the previous subcutaneous injection scenario that's causing the 3+ hour runtime:

1. Tiny Initial Particle Radius (r0) Creates Stiff Dynamics #

Your particle radius r0 = 9.75e-8 m (that's 0.0975 micrometers) is orders of magnitude smaller than typical particles used in drug dissolution models. This tiny value makes the coefficient (D*Cs)/(rho*r0^2) in your first ODE term explode—since r0 is squared in the denominator, the coefficient becomes enormous.

This leads to stiffness in your ODE system: one variable changes extremely rapidly compared to the other. The ode113 solver is designed for non-stiff problems; it will grind to a halt trying to take tiny steps to capture the fast-changing dynamics, while also needing to cover your long simulation window. The previous subcutaneous model almost certainly used larger particles, avoiding this stiffness issue entirely.

2. Extremely Long Simulation Time Span #

You're simulating 200 hours (720,000 seconds) of dissolution—this is an unusually long time frame for an ODE simulation, even for non-stiff systems. Combine this with the stiffness from the tiny particles, and ode113 has to take millions of tiny steps to both resolve the fast initial changes and trudge through the full 200 hours. The prior model likely used a much shorter simulation window, cutting down on total computation time.

3. Missing Absorption Term Alters System Feedback #

You removed the Af*y(2) absorption term from the original model. In the subcutaneous scenario, this term acted as a "sink" for the dissolved drug concentration y(2), keeping the system's dynamics more balanced. Without it, y(2) will continuously rise, creating a stronger nonlinear interaction with y(1) (the normalized particle radius). This increased nonlinearity forces the solver to take more steps to maintain accuracy, adding to the runtime.

4. Astronomically Large Particle Count (N) Amplifies Stiffness #

Calculating N = W/((4/3)*pi*r0^3*rho) gives you a particle count in the 10^15 range—that's quadrillions of particles. This makes the coefficient in your second ODE (D*4*pi*N*r0/V) massive, meaning y(2) changes extremely quickly compared to y(1). This widens the gap between the system's fast and slow dynamics, making the stiffness even worse for ode113 to handle.

Quick Fixes to Speed Up Your Simulation #

• Switch to a stiff solver: Use ode15s or ode23s instead of ode113—these solvers are optimized for stiff systems and will cut down runtime drastically.
• Shorten your simulation window: Check if you really need 200 hours. You can stop the simulation early once y(1) (normalized radius) approaches 0 (i.e., all particles are dissolved).
• Scale your variables: Consider re-scaling your equations to balance the coefficients, reducing the numerical disparity between your ODE terms.

内容的提问来源于stack exchange，提问作者Engineer_1331