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 DDElingen
%DDELINGEN Solve a DDE model. 
%   DDELINGEN compares numerical solution vs. analytical solution of
%   a DDE model with a linear history function.   

tau = 5;
p.a = 1;
p.lag = tau;
p.k = 0.1;
tspan = [0 10];
options = ddeset('RelTol',1e-6, 'AbsTol', 1e-6);
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 Linear History Function');
print('-r600', '-dtiff', 'ddelingen')  


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

% Derivatives vector of DDE system

xlag = Z(1,1);
dxdt = -p.k*xlag;


function xhist = history(t,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.a*t;


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.k*p.a*(-t.^2/2 + p.lag*t) .* (t<p.lag) + ...
    (p.k^2*p.a*(1/6*(t-p.lag).^3 - 1/2*p.lag*(t-p.lag).^2) ...
    + 1/2*p.k*p.a*p.lag^2) .* (p.lag <= t & t <= 2*p.lag);  


Contact us