I am trying to run this code for parameter estimation, but unable to fix the error. please could you help me to fix the error?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Tem_VL
function Tem_VL
% This function fits the first set of BSfludata to and SIR model
clear all
close all
clc
% loading VL data
load VL.txt
% specifying higher precision
format short
tdata = VL(:,1); % define array with t-coordinates of the data
qdata = VL(:,2); % define array with y-coordinates of the data
tforward = 3:0.01:14; % t mesh for the solution of the differential equation
tmeasure = [1:100:1101]; % selects the points in the solution
% corresponding to the t values of tdata
y0=[8381 2000 200 814 800 8029 700 100 14291400 500 150 5000 950]
% initial values of parameters to be fitted
b1 = 0.56;
b2 = 0.75;
b3 = 0.70;
b4 = 0.56;
c1 = 4;
c2 = 8;
c3 = 0.142;
c4 = 0.333;
d1 = 0.018;
d2 = 0.018;
e = 0.33;
f = 0.33;
h1 = 0.58;
h2 = 0.62;
n1 = 6;
n2 = 6;
n3 = 7.5;
n4 = 8.33;
t1 = 0.5294;
t2 = 0.5294;
u = 0.0000367;
v = 0.00714;
m = 0.16;
w1 = 0.104;
w2 = 0.104;
x = 0.000228;
r = 0.16;
z = 0.8;
p=0.95;
j=0.65;
function dy = model_1(t,y,k) % DE
% Assignes the parameters in the DE the current value of the parameters
b1 = k(1);
b2 = k(2);
b3 = k(3);
b4 = k(4);
c1 = k(5);
c2 = k(6);
c3 = k(7);
c4 = k(8);
d1 = k(9);
d2 = k(10);
e = k(11);
f = k(12);
h1 = k(13);
h2 = k(14);
n1 = k(15);
n2 = k(16);
n3 = k(17);
n4 = k(18);
t1 = k(19);
t2 = k(20);
u = k(21);
v = k(22);
m = k(23);
w1 = k(24);
w2 = k(25);
x = k(26);
r = k(27);
z = k(28);
p = k(29);
j = k(30);
% assigns zeros to dy
dy = zeros(14,1);
%NH=y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10)
dy(1)=n1-(u+(b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)...
+y(9)+y(10)))*y(1);% RHS of first equation
dy(2)=((c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)...
+y(10)))*y(1)-(u+r)*y(2);% RHS of second equation
dy(3)=p*r*y(2)-(e+w1+u)*y(3);
dy(4)=r*(1-p)*y(2)+e*y(3)-(t1+d1+u)*y(4);
dy(5)=w1*y(3)+t1*y(4)-u*y(5);
dy(6)=n2-(u+(c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)...
+y(9)+y(10)))*y(6);
dy(7)=((c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)...
+y(10)))*y(6)-(u+m)*y(7);
dy(8)=m*j*y(7)-(f+w2+u)*y(8);
dy(9)=m*(1-j)*y(7)+f*y(8)-(t2+d2+u)*y(9);
dy(10)=w2*y(8)+t2*y(9)-u*y(10);
dy(11)=n3-(v+c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)...
+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)+y(14)))*y(11);
dy(12)=(c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)+y(3)+y(4)...
+y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)...
+y(14)))*y(11)-v*y(12);
dy(13)=n4-(x+(c2*b4*y(12))/(y(13)+y(14)))*y(13);
dy(14)=((c2*b4*y(12))/(y(13)+y(14)))*y(13)-x*y(14);
end
%computing the error in the data
function error_in_data = moder(k)
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
% solves the DE; output is written in T and Y
% assignts the y-coordinates of the solution at the t-coordinates of tdata
q = Y(tmeasure(:),2);
%computes SSE
error_in_data = sum((q - qdata).^2)
end
% main routine; assigns initial values of parameters
k = [b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 e f h1 h2 n1 n2 n3 n4 ...
t1 t2 u v m w1 w2 x r z p j];
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
% solves the DE with the initial values of the parameters
yint = Y(tmeasure(:),2);
% assigns the y-coordinates of the solution at tdata to yint
figure(1)
subplot(1,2,1);
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'g-'); % plotting of solution and data with initial
% guesses for the parameters
xlabel('time in days');
ylabel('Number of cases');
axis([3 14 0 350]);
[k,fval] = fminsearch(@moder,k); % minimization routine; assigns the new
% values of parameters to k and the SSE
% to fval
disp(k);
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
% solving the DE with the final values of the
% parameters
yint = Y(tmeasure(:),2); % computing the y-coordinates corresponding to the
% tdata
subplot(1,2,2)
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'c-');
xlabel('time in days'); % plotting final fit
ylabel('Number of cases');
axis([3 20 0 350]);
end
4 comentarios
Respuestas (1)
Walter Roberson
el 12 de En. de 2023
y0=[8381 2000 200 814 800 8029 700 100 14291400 500 150 5000 950]
That is 13 entries, one of which is much larger value than the rest. We might wonder, for example, whether the data is intended to be
y0=[8381 2000 200 814 800 8029 700 100 1429 1400 500 150 5000 950]
??
dy = zeros(14,1);
Your code is obviously expecting to work with 14 state variables, but you are only passing in 13 state values.
3 comentarios
Walter Roberson
el 12 de En. de 2023
tdata = VL(:,1); % define array with t-coordinates of the data
That is whatever size is read from the file.
tforward = 3:0.01:14
That is probably 1101 elements.
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
You pass those 1101 elements as the times to return for ode23s. Those will be returned into T (unless the ode23s has to end early)
tmeasure = [1:100:1101]; % selects the points in the solution
That is 11 points.
yint = Y(tmeasure(:),2); % computing the y-coordinates corresponding to the
That is going to be 11 points.
plot(tdata,yint,'g-'); % plotting of solution and data with initial
tdata is whatever size was read from the file. yint is 11 points selected out of the 1101 generated for the ode. Those are probably not going to be the same size.
The values in tdata could potentially not intersect 3:0.01:14 at all, or might overlap, or might be fully contained inside (as outside observers we do not know.) You could consider interpolating the values in tdata against the values in tforward to find the closing matching indices and extract the corresponding Y(:,2) values, and plotting tdata against those values.
Or in the case where you are certain that tdata (from the file) is entirely contained within 3:0.01:14 you could create tforward as being [3 tdata] and then throwing away the first row of result (corresponding to the time 3), and the remanining rows would correspond to data. Unless, that is, tdata happens to start with 3 itself. Or unless it happens to start before 3...
Ver también
Categorías
Más información sobre Least Squares en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!