Modellering af fysiologiskesystemer, Home Assginment II

ex1.m

% Home assignment II, exercise 1: Lyapunov exponents
% 
% Emil Munthe, 7. April 2014
%
function [] = ex1(a)
%% Init
if nargin < 1, a = 0.2; end
t = linspace(0, 50, 1500);
yzero = [0.1 0 0];

%% Solve
y = roesler(a,t,yzero);

%% Plot x,y plane
figure(1)
plot(y(1,:),y(2,:))
print(gcf,'-depsc2',['phaseplot' '_a=' num2str(a*100) '.eps'])

%% Liapunov coeefficients
yzero(1) = yzero(1) + 1e-5;
yl = roesler(a,t,yzero);
delta = distance(y,yl);
[r,m,b] = regression(t,log(delta));
figure(2)
plot(t, log(delta),...
     t, m*t+b)
print(gcf,'-depsc2',['deltaplot' '_a=' num2str(a*100) '.eps'])
fprintf('Slope of ln || delta(t) || with a = %f is %f \n ',a,m)

end

function delta = distance(y1,y2)
    delta = sqrt(sum((y2-y1).^2));
end

function [y] = roesler(a,t,yzero)
b = 0.2;
c = 5.7;

if nargin < 2,t = linspace(0, 50, 1500); end
tspan = [0 200];
if nargin < 3, yzero = [0.1
                        0
                        0];
end

% Integration
sol = ode45(@roesler_diff, tspan, yzero);
% Interpolation
y = deval(sol,t,[1 2 3]);

function xdiff = roesler_diff(t,x)
    xdiff = [-x(2)-x(3)
             x(1) + a*x(2)
             b + x(3)*(x(1)-c)];
end

end

ex2.m

% Home assignment II, exercise 2: Medicine concentration
% 
% Emil Munthe, 8. April 2014
% 
%% Init
w = 80;
VA = 0.95 * w;
alpha = 0.28;

%% a)
m = 500e-6;
c0 = m / VA;
fprintf('a) Initial plasma concentration is %e\n',c0)

h = 8;
p = 6;
t = 0:1/60:h;
c = 0;
for i = 1:p
    c_tmp = (c(end)+c0)*exp(-alpha*t);
    c = [c(1:end-1) c_tmp];
end
t = linspace(0,h*p,length(c));
figure(1)
plot(t,c)

%% b)
m = 125e-6;
c0 = m / VA;
fprintf('b) Initial plasma concentration is %e\n',c0)

h = 2;
p = 24;
t = 0:1/60:h;
c = 0;
for i = 1:p
    c_tmp = (c(end)+c0)*exp(-alpha*t);
    c = [c(1:end-1) c_tmp];
end
t = linspace(0,h*p,length(c));
figure(1)
hold on
plot(t,c,'r')
hold off

%% c)
%Numbers from b), measured on the plot
c_max = 3.818e-6;
c_min = 2.201e-6;

m = 285e-6;
c0 = m / VA;
md = 125e-6;
cd = md /VA;
h = 2;
p = 24;
t = 0:1/60:h;
% Initial dosis:
c = c0*exp(-alpha*t);
% Cyclic dodis, cd:
for i = 1:p
    c_tmp = (c(end)+cd)*exp(-alpha*t);
    c = [c(1:end-1) c_tmp];
end
t = linspace(0,h*p,length(c));
fprintf('Min of concentration in b) %e and in c) %e\n',c_min, min(c))
fprintf('Max of concentration in b) %e and in c) %e\n',c_max, max(c))

figure(1)
hold on
plot(t,c,'g')
hold off

%% d)
t = [10 20 30 45 60];
c = [29.6 17.8 12.6 10.4 6.0];

m = 30e-6;
[r,a,b] = regression(t,log(c));
VA = (m/exp(b))*1e6;

figure(2)
plot(t,log(c))

fprintf('Rats have a apparent volume of distribution of %f l/kg and elimination rate constant of %f /h \n',VA,-a)