Code covered by the BSD License  

Highlights from

image thumbnail




15 Mar 2013 (Updated )

MATLAB routines for the book: "Development of Innovative Drugs via Modeling with MATLAB".

function DDEconst
%DDECONST Solve a DDE model. 
%   DDECONST compares numerical solution vs. analytical solution of
%   a DDE model with a constant history.   

tau = 5;
p.x0 = 10;
p.lag = tau;
p.k = 0.1;
tspan = [0 10];
options = ddeset('RelTol',1e-9, 'AbsTol', 1e-9);
sol  = dde23(@derivatives, p.lag, @history, tspan, options, p);
% numerical solution
t  = sol.x;
x  = sol.y; 
% ax = axes;
plot(t,x,'--b','LineWidth',2); hold on
% numerical solution
x = analyticsol(t,p);
legend('numerical', 'analytical')
title('DDE Solution with Constant History Function');
print('-r600', '-dtiff', 'ddeconst')  


function dxdt = derivatives(~, ~, Z, p)
%DERIVATIVES Compute the right-hand side of the DDE.
%   DYDT = DERIVATIVES(T, Y, Z, P) calculates |DYDT|, he right-hand
%   side of the DDE model, at points defined by the vector |Y|, |T|,
%   and with parameters |P|. Matrix |Z| corresponds to all possible
%   delays and contains solutions |Y| at the delayed argument. Thus,
%   |Z| has the number of rows equal to the length of |Y|, and the
%   number of column equal to the number of delays. For example:
%   Z(3,1) means y(3) at time |T|-lag; lag is a delay. 

% Derivatives vector of the DDE system
xlag = Z(1,1);
dxdt = -p.k*xlag;


function xhist = history(~,p)
%HISTORY Provide the history of the solution 
%   XHIST = HISTORY(T, P) specifies the time profiles of dependent
%   variables at times |T| prior to the initial point. This function
%   is a form of 'initial condition' for the DDE. |P| stands for
%   additional parameters which can be passed to the function.

xhist = p.x0;


function x = analyticsol(t,p)
%ANALYTICSOL Analytical solution of the DDE
%   X = ANALYTICSOL(T,P)computes the dependent variable |X|, defined
%   at time |T|, and based on the analytical solution with 
%   parameters |P|.

x = p.x0*(1-p.k*t) .* (t <= p.lag) +                ... 
    p.x0*(1 - p.k*t + (p.k^2)/2*(t-p.lag).^2) .*    ...
    (p.lag < t & t <= 2*p.lag);  


Contact us