Figure A.3:

Solution to first-order differential equation dc_A/dt=-kc_A, and sensitivities S_1=d c_A/d k and S_2=d c_A/d c_{A0}.

Code for Figure A.3

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

%
% Illustrate sensitivity calculation in simple first-order reaction
%
% dca/dt = -k ca
% ca(0)  = ca0
%
% theta = [k; ca0]
%
% compute solution to model
% and sensitivities
%
%  ca(t)
%  S1 = d ca / d k
%  S2 = d ca / d ca0
%
%  ca(t) = ca0 exp(-kt)
%  S1 = -t ca0 exp(-kt)
%  S2 = exp(-kt)
%
% jbr,  10/5/01
%
% converted from ddasac to cvodes, 3/31/07, jbr
%
% converted to use paresto, 5/1/2018, jbr
%

ca0 = 2;
k   = 1;
tfinal = 5;
nts    = 100;
tout = linspace(0,tfinal,nts);

model = struct;
model.transcription = 'shooting';
model.x = {'ca'};
model.p = {'k'};
model.d = {'m_ca'};
model.tout = tout;

model.ode = @(t, y, p) {-p.k*y.ca};
% fit the perfect measurement data
model.lsq = @(t, y, p) {y.ca-y.m_ca};

pe = paresto(model);

% initialize state and parameters

x0 = ca0;
thetaac = k;

y_ac = pe.simulate(0, x0, thetaac);

theta0 = struct();
theta0.k = k;
theta0.ca = ca0;

lb.k = theta0.k;
lb.ca = theta0.ca;

ub.k = theta0.k;
ub.ca = theta0.ca;

[est, y, p] = pe.optimize(y_ac, theta0, lb, ub);

table = [tout(:), y.ca(:), est.dca_dk(:), est.dca_dca0(:)];
save -ascii Sfirstorder.dat table;
% if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (tout, [y.ca', est.dca_dk, est.dca_dca0]);
% % TITLE
% end % PLOTTING