estimate the relative sensitivity coefficients
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/911130/image.png)
% TODO - Write the function declaration. Name the function SIRmodel
function dYdt = Lab3(t,Y,p)
% TODO - Extract S1, S2
s1 = Y(1);
s2 = Y(2);
% TODO - Define the constants/parameters in the ODE system
V0 = p(1);
k1 = p(2);
V2 = p(3);
KM = p(4);
ds1dt = V0 - k1*s1;
ds2dt = k1*s1 - ((V2*s2)/(KM+s2));
% TODO - Create output column vector dYdt
dYdt = [ds1dt; ds2dt];
end
clear all
% Close all
% Timespan
tRange = [0 100];
% Initial condition
Y0 = [0;0];
% Parameters in ODE system
V0 = 1; % mM/min
k1 = 25; % 1/min
V2 = 2; % mM/min
KM = 5; % mM
% Want to choose a parameter p and a value between 1% and 5% so set vector
% to be P*[1;1.03]; so this give back the parameter the parameter with a 3% increase
delp = 0.01;
% For changing the param values in the loop
p_delp = V0*[1;1+delp];
% Run the model for the two values of the parameter p and p+delta_p
for i=1:length(p_delp)
% Choose p and p + delta_p estimating the sensitivity
V0 = p_delp(i);
% Set the parameters to send to solver
p = [V0,k1,V2,KM];
% Call the ODE solver ode15s instead of ode45
% to send parameters to the ode solver, use the following command:
[tSol,YSol] = ode15s(@(tSol,YSol)Lab3(tSol,YSol,p),tRange,Y0);
s1 = YSol(:,1);
s2 = YSol(:,2);
% Extract the steady state
s1ss(i) = s1(end);
s2ss(i) = s2(end);
end
% Compute dS1ss/dV0 and dS2ss/dV0
ds1ssdV0 = (s1ss(2) - s1ss(1))/delp;
ds2ssdV0 = (s2ss(2) - s2ss(1))/delp;
% Rescale the derivatives:
% new: (V0/s1ss)*ds1ss/dV0
% new: (V0/s2ss)*ds2ss/dV0
% Coefficients
scs1 = (V0./s1ss)*ds1ssdV0
scs2 = (V0./s2ss)*ds2ssdV0
% Plot solutions in time
figure(1)
clf
hold on
plot(tSol,YSol,'LineWidth',2)
xlabel('Time')
ylabel('Concentration')
legend('S_1','S_2')
set(gca,'FontSize',18)
grid on
0 Kommentare
Antworten (1)
Torsten
am 1 Mär. 2022
% TODO - Write the function declaration. Name the function SIRmodel
tRange = [0 100];
% Initial condition
Y0 = [0;0];
% Parameters in ODE system
V0 = 1; % mM/min
k1 = 25; % 1/min
V2 = 2; % mM/min
KM = 5; % mM
% Want to choose a parameter p and a value between 1% and 5% so set vector
% to be P*[1;1.03]; so this give back the parameter the parameter with a 3% increase
delp = 0.01;
% For changing the param values in the loop
p_delp = V0*[1;1+delp];
% Run the model for the two values of the parameter p and p+delta_p
for i=1:length(p_delp)
% Choose p and p + delta_p estimating the sensitivity
V0 = p_delp(i);
% Set the parameters to send to solver
p = [V0,k1,V2,KM];
% Call the ODE solver ode15s instead of ode45
% to send parameters to the ode solver, use the following command:
[tSol,YSol] = ode15s(@(tSol,YSol)Lab3(tSol,YSol,p),tRange,Y0);
% Extract the steady state
s1ss(i) = YSol(end,1);
s2ss(i) = YSol(end,2);
% Save results
t{i} = tSol;
s1{i} = YSol(:,1);
s2{i} = YSol(:,2);
end
% Compute dS1ss/dV0 and dS2ss/dV0
ds1ssdV0 = (s1ss(2) - s1ss(1))/delp;
ds2ssdV0 = (s2ss(2) - s2ss(1))/delp;
% Rescale the derivatives:
% new: (V0/s1ss)*ds1ss/dV0
% new: (V0/s2ss)*ds2ss/dV0
% Coefficients
scs1 = (p_delp(1)/s1ss(1))*ds1ssdV0
scs2 = (p_delp(1)/s2ss(1))*ds2ssdV0
% Plot solutions in time
figure(1)
plot(t{1},[s1{1},s2{1}],'LineWidth',2)
hold on
plot(t{2},[s1{2},s2{2}],'LineWidth',2)
xlabel('Time')
ylabel('Concentration')
legend('S_1','S_2')
set(gca,'FontSize',18)
grid on
end
0 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!