drugdelivery.m

这是一些关于MATLAB的小程序

function DrugDelivery
% [Ref] Problem Solving in Chemical Engineering
% by Michael B. Cutlip and Mordechai Shacham.
% 7.4(a) Controlled Drug Delivery by Dissolution of Pill Coating
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/06/26 $

clear all
clc
global rho SA SD V tau
rho = 1414.7;
SA = 1.0;
SD = 0.4;
V = 1200;
tau = 240;
t0 = 0;
tf = 150;
D10 = 0.5;
D20 = 0.4;
D30 = 0.35;
CAS0 = 0;
CDS0 = 0;
tspan = [t0 tf];
y0 = [D10 D20 D30 CAS0 CDS0];
[t,y] = ode45(@ODEs,tspan,y0)

% 输出结果
plot(t,y(:,1),'-',t,y(:,2),'.-',t,y(:,3),'--')
xlabel('时间 (min)')
ylabel('药丸直径 (cm)')
legend('D_1','D_2','D_3')
title('三颗药丸在溶解过程中的直径变化')
figure
plot(t,y(:,4),'-',t,y(:,5),'--')
xlabel('时间 (min)')
ylabel('胃中的药物浓度 (mg/cm^3)')
legend('C_A_S','C_B_S')
title('胃中药物浓度的动态变化')


% ------------------------------------------------------------------
function dydt = ODEs(t,y)
global rho SA SD V tau
D1 = y(1);
D2 = y(2);
D3 = y(3);
CAS = y(4);
CDS = y(5);
kL1 = 2*0.6/D1;
kL2 = 2*0.6/D2;
kL3 = 2*0.6/D3;
if D1>0.3
    SW1 = 1.0;
else
    SW1 = 0.0;
end
if D2>0.3
    SW2 = 1.0;
else
    SW2 = 0.0;
end
if D3>0.3
    SW3 = 1.0;
else
    SW3 = 0.0;
end

if D1>0.3
    dD1dt = -2*kL1*(SA-CAS)/rho;
elseif D1>1e-6
    dD1dt = -2*kL1*(SD-CDS)/rho;
else
    dD1dt = 0.0;
end

if(D2>0.3)
    dD2dt = -2*kL2*(SA-CAS)/rho;
elseif D2>1e-6
    dD2dt = -2*kL2*(SD-CDS)/rho;
else
    dD2dt = 0.0;
end

if D3>0.3
    dD3dt = -2*kL3*(SA-CAS)/rho;
elseif D3>1e-6
    dD3dt = -2*kL3*(SD-CDS)/rho;
else
    dD3dt = 0.0;
end

dCASdt = (1/V)*(SW1*kL1*(SA-CAS)*pi*D1^2+SW2*kL2*(SA-CAS)*pi*D2^2  ...
        +SW3*kL3*(SA-CAS)*pi*D3^2)-CAS/tau;
dCDSdt = (1/V)*((1-SW1)*kL1*(SD-CDS)*pi*D1^2+(1-SW2)*kL2*   ...
        (SD-CDS)*pi*D2^2+(1-SW3)*kL3*(SD-CDS)*pi*D3^2)-CDS/tau;
dydt = [dD1dt; dD2dt; dD3dt; dCASdt; dCDSdt];