estimation of parameters by using lsqnonlin for system of differential equations
8 次查看(过去 30 天)
显示 更早的评论
1 36456
2 33761
3 38309
4 36638
5 36111
6 32272
7 27107
8 32067
9 26357
10 34666
11 30034
12 30354
13 27336
14 21941
15 26401
16 18164
17 26762
18 26991
19 26834
20 24589
21 19174
22 23886
23 24236
24 23924
25 22350
26 18574
27 20333
28 16072
29 20529
30 21957
31 19046
32 19522
33 19460
34 18938
35 18967
36 18591
37 18381
38 18252
39 18008
40 18105
41 18020
42 17478
43 17191
44 17297
45 16348
46 15877
47 15338
48 15034
49 14678
50 14376
51 14128
52 13959
53 13840
54 13818
55 13714
56 13588
57 13437
58 12905
59 13535
60 13354
61 13094
62 12850
63 12776
64 12529
65 12724
66 11794
67 11602
68 11459
69 11496
70 11515
71 11477
72 11422
73 10986
74 11047
75 11066
76 11054
77 11107
78 11230
79 11276
80 11831
81 12085
82 12331
83 12699
84 12900
85 13198
86 13780
87 14260
88 14638
89 15050
90 15241
91 15495
92 15682
93 15753
94 15785
95 16037
96 16311
97 16745
98 17185
99 17596
100 18371
101 19296
102 20228
103 21147
104 22270
105 23567
106 25138
107 26994
108 29335
109 31629
110 34295
111 37223
112 39537
113 42161
114 44673
115 47444
116 50497
117 53185
118 56213
119 58430
120 59287
121 61958
122 65145
123 68966
124 73304
125 78388
126 84160
127 93028
128 100766
129 107977
130 115974
131 124484
132 133929
133 143114
134 153113
135 163586
136 175723
137 188437
138 203913
139 218937
140 232675
141 248257
142 264852
143 281380
144 297279
145 309865
146 321233
147 330155
148 339945
149 349038
150 356787
151 364870
152 371115
153 373320
154 378505
155 381353
156 386099
157 390029
158 390727
159 392331
160 391812
161 388058
162 383159
163 376017
164 342794
165 354315
166 341022
167 328934
168 319438
169 307822
170 295473
171 283507
172 273656
173 263676
174 255247
175 245652
176 237330
177 228090
178 217638
179 205750
180 194948
181 185028
182 175174
183 142207
184 153274
185 145609
186 137948
187 130692
188 123236
189 117368
190 111601
191 105864
192 100068
193 94942
194 90090
195 85775
196 82090
197 77722
198 73923
199 69721
200 66320
201 63190
202 60615
203 58140
204 56512
205 54658
206 53118
207 51404
208 50151
209 49229
210 48427
211 47754
212 46939
213 46242
214 45588
215 44614
216 43704
217 43269
218 42963
219 42548
220 42080
221 41862
222 41643
223 41286
224 40827
225 40306
226 39743
227 39110
228 38461
229 38431
230 38527
231 38328
232 38587
233 38577
234 37975
235 38173
236 38031
237 38009
238 38209
239 38330
240 38541
241 39942
242 40227
time = data(:,1);
A_data= data(:,2);
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@Kinetics,B0,time,A_data,lb,ub,options );
disp(B)
function Av = Kinetics(B,time);
plot(time,A_data,time,Av)
function A = Kinetics(B, t)
x0 = [1217378052,100,3,2,1,1,1,1];
[T,Av] = ode45(@DifEq, t, x0);
function dA = DifEq(t, x)
N = 1390000000;
pi = 70000;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,1);
end
i got an error;
Error using lsqnonlin (line 190)
Invalid datatype. Options argument must be created with OPTIMOPTIONS.
Error in error_14bmodel (line 19)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@Kinetics,B0,time,A_data,lb,ub,options );
i am fresher of matlab .i requesting you sir please help anyone how to resolve this error
if possible please refer some covid 19 or malaria,dengue hepities b models sir
thank you sir
1 个评论
回答(4 个)
Torsten
2022-4-2
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@Kinetics,B0,time,A_data,lb,ub,options );
This is not the correct call to lsqnonlin, but to lsqcurvefit.
3 个评论
Torsten
2022-4-2
You modify this by calling lsqcurvefit instead of lsqnonlin:
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics,B0,time,A_data,lb,ub,options );
mallela ankamma rao
2022-4-2
Thank you sir I request you sir Please modify this by Lsqnonlin .i didn't know how to modify this.i am struggling alot sir
Star Strider
2022-4-2
I recognise my code.
It should have the correct lsqcurvefit call, so all you need to do is substitute your differnetial equations and data into it.
An updated version of that code is:
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835];
theta0=rand(6,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
fprintf(1,'\tRate Constants:\n')
Rate Constants:
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
Theta(1) = 0.89692
Theta(2) = 0.25742
Theta(3) = 0.31073
Theta(4) = 0.61683
Theta(5) = 0.62818
Theta(6) = 0.00000
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
function C=kinetics(theta,t)
c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dC=dcdt;
end
C=Cv;
end
.
37 个评论
mallela ankamma rao
2022-4-2
Thanks alot for your reply sir Yes i used your code sir. I requesting you sir Please tell me how to apply Lsqnonlin for this code. I tried alot but i couldn't get sir.please help me
Star Strider
2022-4-2
Do not use lsqnonlin. Use lsqcurvefit just as in my code.
If you absolutely must use lsqnonlin, the objective function to use with it is:
fun = @(theta)norm(kinetics(theta,t) - c);
That is the ‘fun’ argument to it in the lsqnonlin documentation.
This references my code. Make approriate changes to use it with your data and your ‘kinetics’ function.
.
Torsten
2022-4-2
I think the function to be supplied for lsqnonlin is
fun =@(theta) kinetics(theta,t) -c
to make the objective equal to the one for lsqcurvefit.
mallela ankamma rao
2022-4-2
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics,B0,time,A_data,lb,ub,options );
Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.
<stopping criteria details>
1.0000
0.0000
0.0000
0.0000
0.0000
1.0000
1.0000
1.0000
0.0000
0.0000
for lsqcurvefit i got very poor parameters sir
so many mathematical modeling on covid 19 papers prefers only lsqnonlin so that i asked how to change this sir
i request you sir please send one line code like this [B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics,B0,time,A_data,lb,ub,options ) for my data
for lsqnonlin
Star Strider
2022-4-3
Try this —
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(theta) kinetics(theta,t) - c, B0,time,A_data,lb,ub,options );
@Torsten — Correct. I rarely use lsqnonlin (and have not in years). The objective function needs to return an array of differences.
.
mallela ankamma rao
2022-4-4
respected sir
i have done by modfying these but i got same values i requesting you sir please help me in this
i have tried many times but i didnot get correct values
i requesting you sir please spend 10 minutes time for me by running this and give me final code sir because i have to apply all my data to this lsqnonlin method.only
thanks alot sir
A_data= [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
tspan = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ];
N = 13900000;
e0 = 121737;
a0 = 150;
i0 = 100;
q0 = 30;
h0 = 10;
r0 = 10;
d0 = 0;
s0 = N-e0-a0-i0-q0-h0-r0-d0;
x0 = [s0; e0; a0; i0; q0; h0; r0; d0];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(B,A_data,tspan,x0),B0,lb,ub,options );
disp(B)
function A = immanuel(B,tspan,x0)
[T,Av] = ode45(@(t,x)DifEq(t, x),tspan, x0);
function dA = DifEq(t, x)
N = 1390000000;
pi = 70000;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,2);
end
ERROR
Error using lsqnonlin_3data1>immanuel
Too many input arguments.
Error in lsqnonlin_3data1>@(B)immanuel(B,A_data,tspan,x0) (line 28)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(B,A_data,tspan,x0),B0,lb,ub,options );
Error in lsqnonlin (line 218)
initVals.F = feval(funfcn{3},xCurrent,varargin{:});
Error in lsqnonlin_3data1 (line 28)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(B,A_data,tspan,x0),B0,lb,ub,options );
Caused by:
Failure in initial objective function evaluation. LSQNONLIN cannot continue.
----------------------------------------------------------------------------------------------------------------
clc;
A_data= [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
tspan = [1 2 3 4 5 7 9 11 13 15 17 21 25 30 40 50 60];
N = 1390000000;
e0 = 121737;
a0 = 150;
i0 = 100;
q0 = 30;
h0 = 10;
r0 = 10;
d0 = 0;
s0 = N-e0-a0-i0-q0-h0-r0-d0;
x0 = [s0; e0; a0; i0; q0; h0; r0; d0];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(A_data,tspan,x0),B0,lb,ub,options );
disp(B)
function A = immanuel(B,tspan,x0)
[T,Av] = ode45(@(t,x)DifEq(t, x),tspan, x0);
function dA = DifEq(t, x)
N = 1390000000;
pi = 70000;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,2);
end
i got values similar values
Initial point is a local minimum.
Optimization completed because the size of the gradient at the initial point
is less than the value of the optimality tolerance.
<stopping criteria details>
0.5000
0.0040
0.1000
0.1300
0.0700
0.1000
0.0700
0.0300
0.0001
0.0002
Torsten
2022-4-4
编辑:Torsten
2022-4-5
The list of inputs for your function "immanuel" (Kant ?) must coincide in the call to lsqnonlin and in the function itself.
In your first attempt
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(B,A_data,tspan,x0),B0,lb,ub,options );
disp(B)
function A = immanuel(B,tspan,x0)
you tell MATLAB that the list of inputs to "immanuel" is
@(B)immanuel(B,A_data,tspan,x0)
but the function itself is defined as
function A = immanuel(B,tspan,x0)
Thus the argument "A_data" in the second position is missing.
In your second attempt where the initial values are returned, you make a similar mistake.
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(A_data,tspan,x0),B0,lb,ub,options );
disp(B)
function A = immanuel(B,tspan,x0)
Do you see the difference in the input lists
@(B)immanuel(A_data,tspan,x0)
and
immanuel(B,tspan,x0)
?
Furthermore, since you use lsqnonlin and not lsqcurvefit,
A = A_data - Av(:,2)
and not
A = Av(:,2)
is the return parameter from "immanuel" if your experimental data are to be compared with the second solution variable of the ODE integrator.
mallela ankamma rao
2022-4-5
its worked sir.thanks alot Torston sir
and thanks alot star strider sir
mallela ankamma rao
2022-4-5
respected sir i want to fix the same data for using lsqcurvefit but i could not get values it shows error sir
sir please guide me in this method also where i did wrong and please modify that sir
thank you sir
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions('lsqcurvefit', 'Display', 'iter','Algorithm', 'levenberg-marquardt','UseParallel', false,'StepTolerance', 1e-6);
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
disp(B)
Av = Kinetics(B,t);
plot(t,x,t,Av)
function A = Kinetics(B, t)
x0 =[1217378,12173,300,200,100,100,100,90];
[T,Av] = ode45(@DifEq, t, x0);
function dA = DifEq(t, x)
N = 139000000;
pi = 1500;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av;
end
ERROR
Error using lsqcurvefit (line 286)
Function value and YDATA sizes are not equal.
Error in error_14bmodel (line 17)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
Star Strider
2022-4-5
The problem is that as written, ‘A’ returns a matrix. You are only fitting one column with the output, so choose whatever column that is, for example:
A = Av(:,1); % Return Only One Column To Fit A Column Vector
Thsi arbitrarily chooses the first column, however that may not be correct. Choose the correct column and the code should run with out error and give a reasonable fit (if everything else is coded correctly).
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
x = x(:);
t = t(:);
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions('lsqcurvefit', 'Display', 'iter','Algorithm', 'levenberg-marquardt','UseParallel', false,'StepTolerance', 1e-6);
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
disp(B)
Av = Kinetics(B,t);
plot(t,x,t,Av)
function A = Kinetics(B, t)
x0 =[1217378,12173,300,200,100,100,100,90];
[T,Av] = ode45(@DifEq, t, x0);
function dA = DifEq(t, x)
N = 139000000;
pi = 1500;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,1); % Return Only One Column To Fit A Column Vector
end
.
mallela ankamma rao
2022-4-5
sir i modifiied as you said and applied but i got same error sir
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions('lsqcurvefit', 'Display', 'iter','Algorithm', 'levenberg-marquardt','UseParallel', false,'StepTolerance', 1e-6);
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
disp(B)
Av = Kinetics(B,t);
plot(t,x,t,Av)
function A = Kinetics(B, t)
x0 =[1217378,12173,300,200,100,100,100,90];
[T,Av] = ode45(@DifEq, t, x0);
function dA = DifEq(t, x)
N = 139000000;
pi = 1500;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,2);
end
ERROR
Error using lsqcurvefit (line 286)
Function value and YDATA sizes are not equal.
Error in error_14bmodel (line 18)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
sir if possible please run this code and please change anything if necessary
thank you sir
mallela ankamma rao
2022-4-5
this also same error sir
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions('lsqcurvefit', 'Display', 'iter','Algorithm', 'levenberg-marquardt','UseParallel', false,'StepTolerance', 1e-6);
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
disp(B)
Av = Kinetics(B,t);
plot(t,x,t,Av)
function A = Kinetics(B, t)
x0 =[1217378,12173,300,200,100,100,100,90];
[T,Av] = ode45(@DifEq, t, x0);
function dA = DifEq(t, x)
N = 139000000;
pi = 1500;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,1);
end
Error using lsqcurvefit (line 286)
Function value and YDATA sizes are not equal.
Error in error_14bmodel (line 18)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@Kinetics, B0,t,x, zeros(size(B0)), [], options);
Star Strider
2022-4-5
That is likely caused by the fact that the differential equation integrators return column vectors (or column-oriented matrices).
Force ‘x’ and ‘t’ to be column vectors.
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
x = x(:);
t = t(:);
That will then work.
Also, Torsten believes that:
A = Av(:,2);
is correct. Change the column reference to whatever works with your code to regress against the correct column.
.
Torsten
2022-4-5
Testing your code shows that changing
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
to
x = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346].';
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].';
removes the problem.
mallela ankamma rao
2022-4-6
i got it . Thank you very much star strider sir
Thank you very much torsten sir.
mallela ankamma rao
2022-4-6
good morning sir
with your guidence i learned how to fit parameters by using fmin,lsqnonlin,lsqcurvefit.
thanks alot sir
sir i have seen in your previous codes , ga function also give good fit for parameters
sir please guide me how to apply that ga function to my code.
thank you sir
Star Strider
2022-4-6
A related example —
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835];
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 10;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-8, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
function C=kinetics(theta,t)
c0 = theta(7:10);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dC=dcdt;
end
C=Cv;
end
Its running time exceeds the 55 seconds that the online Run feature enforces, so copy it as is, save it to a new .m file on your MATLAB search path, and run it to see how it works. When I ran it a few minutes ago, it required about 11 minutes to converge (AMD Ryzen 9 3900 processor).
.
mallela ankamma rao
2022-4-6
sir i have tried but it shows error sir. i have one doubt
how to take inital conditions of parameters
theta0 = [0.5,0.004,0.1,0.13,0.07,0.1,0.07,0.03,0.0001,0.0002 ]; in
[p,fval,exitflag,output] = ga(@(p)ftns(p,t),Parameter,[],[],[],[],zeros(Parameter,1),[ones(10,1); Inf(8,1)],[],[],opts);
please guide me sir by changing this code
thank you sir
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
c = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
theta = [0.5,0.004,0.1,0.13,0.07,0.1,0.07,0.03,0.0001,0.0002 ];
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parameter = 18;
opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,10)*1E-3 ones(1,8)*10],'MaxGenerations',2E3);%,'PlotFcn','gaplotbestf'); %how can I interpret the InitialPopulationMatrix?
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n',t0)
[p,fval,exitflag,output] = ga(@(p)ftns(p,t),Parameter,[],[],[],[],zeros(Parameter,1),[ones(10,1); Inf(8,1)],[],[],opts);
disp(p)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0);
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
function C=kinetics(theta,t)
c0 = [1217,100,3,2,1,1,1,1];
[T,CV]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
N = 13900;
pi = 700;
zeta = 0.1;
eta = 0.2;
rita = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = theta(1);
alpha =theta(2);
gamma = theta(3);
upsilon =theta(4);
epsilon = theta(5);
lamda = theta(6);
sigma = theta(7);
kappa = theta(8);
nu = theta(9);
xi = theta(10);
dcdt = zeros(8,1);
dcdt(1) = pi -beta*(zeta*c(3)+eta*c(4)+rita*c(5)+iota*c(6))*(c(1)/N) -mu*c(1);
dcdt(2) = beta*(zeta*c(3)+eta*c(4)+rita*c(5)+iota*c(6))*(c(1)/N) -(delta+mu)*c(2);
dcdt(3) = rho*delta*c(2)-(lamda+gamma+nu+mu)*c(3);
dcdt(4) = (1-rho)*delta*c(2)-(sigma+kappa+mu)*c(4);
dcdt(5) = lamda*c(3) + sigma*c(4)-(alpha+upsilon+mu)*c(5);
dcdt(6) = alpha*c(6) + kappa*c(4)- (epsilon+xi+mu)*c(6);
dcdt(7) = gamma*c(3) + upsilon*c(5) + epsilon*c(6);
dcdt(8) = nu*c(3) + xi*c(6);
dC=dcdt;
end
C=CV;
end
Error using error_sum_ga1>@(theta)norm(c-kinetics(theta,t))
Too many input arguments.
Error in error_sum_ga1>@(p)ftns(p,t) (line 13)
[p,fval,exitflag,output] = ga(@(p)ftns(p,t),Parameter,[],[],[],[],zeros(Parameter,1),[ones(10,1); Inf(8,1)],[],[],opts);
Error in createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)
fcn_handle = @(x) fcn(x,FcnArgs{:});
Error in makeState (line 52)
firstMemberScore = FitnessFcn(state.Population(initScoreProvided+1,:));
Error in galincon (line 22)
state = makeState(GenomeLength,FitnessFcn,Iterate,output.problemtype,options);
Error in ga (line 402)
[x,fval,exitFlag,output,population,scores] = galincon(FitnessFcn,nvars, ...
Error in error_sum_ga1 (line 13)
[p,fval,exitflag,output] = ga(@(p)ftns(p,t),Parameter,[],[],[],[],zeros(Parameter,1),[ones(10,1); Inf(8,1)],[],[],opts);
Caused by:
Failure in initial user-supplied fitness function evaluation. GA cannot continue.
Star Strider
2022-4-6
This seems to run without error —
t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
c = [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = t(:);
c = c(:);
% theta = [0.5,0.004,0.1,0.13,0.07,0.1,0.07,0.03,0.0001,0.0002 ];
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parameter = 18;
opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,10)*1E-3 ones(1,8)*10],'MaxGenerations',1E4, 'FunctionTolerance',1E-10);%, 'PlotFcn','gaplotbestf'); %how can I interpret the InitialPopulationMatrix?
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n',t0)
Start Time: 2022-04-06 15:57:39.3124
[p,fval,exitflag,output] = ga(ftns,Parameter,[],[],[],[],zeros(Parameter,1),[ones(10,1); Inf(8,1)],[],[],opts);
Optimization terminated: maximum number of generations exceeded.
disp(p)
0.9976 0.4964 0.0033 0.0211 0.2662 0.3121 0.5904 0.6113 0.0315 0.2361 138.0790 424.3110 183.2415 104.5209 126.3709 52.9959 24.4418 34.8537
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
Stop Time: 2022-04-06 15:58:10.1306
GA_Time = etime(t1,t0);
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
Elapsed Time: 3.081816300000000E+01 00:00:30.818
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t = %d\n\n',fval,output.generations)
Fitness value at convergence = 501.5223
Generations = 2000
fprintf(1,'\tRate Constants:\n')
Rate Constants:
for k1 = 1:length(p)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, p(k1))
end
Theta( 1) = 0.99756
Theta( 2) = 0.49640
Theta( 3) = 0.00331
Theta( 4) = 0.02109
Theta( 5) = 0.26618
Theta( 6) = 0.31208
Theta( 7) = 0.59045
Theta( 8) = 0.61127
Theta( 9) = 0.03146
Theta(10) = 0.23611
Theta(11) = 138.07903
Theta(12) = 424.31100
Theta(13) = 183.24153
Theta(14) = 104.52093
Theta(15) = 126.37089
Theta(16) = 52.99594
Theta(17) = 24.44183
Theta(18) = 34.85370
tv = linspace(min(t), max(t));
Cfit = kinetics(p, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
function C=kinetics(theta,t)
% c0 = [1217,100,3,2,1,1,1,1];
c0 = theta(11:18);
[T,CV]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
N = 13900;
pi = 700;
zeta = 0.1;
eta = 0.2;
rita = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = theta(1);
alpha =theta(2);
gamma = theta(3);
upsilon =theta(4);
epsilon = theta(5);
lamda = theta(6);
sigma = theta(7);
kappa = theta(8);
nu = theta(9);
xi = theta(10);
dcdt = zeros(8,1);
dcdt(1) = pi -beta*(zeta*c(3)+eta*c(4)+rita*c(5)+iota*c(6))*(c(1)/N) -mu*c(1);
dcdt(2) = beta*(zeta*c(3)+eta*c(4)+rita*c(5)+iota*c(6))*(c(1)/N) -(delta+mu)*c(2);
dcdt(3) = rho*delta*c(2)-(lamda+gamma+nu+mu)*c(3);
dcdt(4) = (1-rho)*delta*c(2)-(sigma+kappa+mu)*c(4);
dcdt(5) = lamda*c(3) + sigma*c(4)-(alpha+upsilon+mu)*c(5);
dcdt(6) = alpha*c(6) + kappa*c(4)- (epsilon+xi+mu)*c(6);
dcdt(7) = gamma*c(3) + upsilon*c(5) + epsilon*c(6);
dcdt(8) = nu*c(3) + xi*c(6);
dC=dcdt;
end
C=CV(:,2);
end
Giving it more iterations and a larger 'InitialPopulationMatrix' that I did here (I originally limited it to 50 rows in order to work with the onlilne Run feature) will llikely provide a better fit to the data.
The first 10 values of ‘theta’ are the parameters, and the last 8 are the initial conditions for the differential equations.
.
mallela ankamma rao
2022-4-7
good morning sir
thank you very much sir for spending your valuable time for me.
it runs good sir.
i have one doubt sir in my paper i have taken the initial conditions
theta = [0.5,0.004,0.1,0.13,0.07,0.1,0.07,0.03,0.0001,0.0002 ] and c0 = [12170,100,30,20,10,10,1,1];
By using lsqnonlin and lscurvfit and fmin we have taken above both values theta and co but in above code it did not take bothe above values.
was it not affect my problem parameters and graph?
just i want to know the reason to explain my guide if he asked the same question sir.
thank you sir
Star Strider
2022-4-7
My code considers the initial conditions, as well as the kinetic constants, as parameters to be estimated. I would use the parameters and initial conditions my code estimates.
It may be necessary to run the genetic algorithm several times (probably at least 10) in order to get the best parameter estimates. Choose the set with the lowest fitness value at convergence.
If they are different from the publised data, they may be correct and the published data may not be. It would be necessary to know how the authors estimated their parameters in order to be certain the published parameters are correct. You can always run the published parameters and initial conditions with your integrated differential equations to see how well they fit the data. I encourage you to do that.
mallela ankamma rao
2022-4-7
ok sir i understood but it takes long time sir .it has been running since 2 hours but not over still.
remaining codes using lsqnonlin and lscurvefit the result came with in 2 minutes.
was it my internet problem or this model takes long time sir?
Star Strider
2022-4-7
The ga function can take a while to converge. It will be necessary to be certain that you are returning the correct column of ‘Cv’ as ‘C’. It converged relatively rapidly when I ran it.
Use this options structure:
opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,10)*1E-3 ones(1,8)*10],'MaxGenerations',1E4, 'FunctionTolerance',1E-10, 'PlotFcn','gaplotbestf');
That will plot the progress of the optimisation so you can see how it is working.
Also, defining:
PopSz = 250;
will speed up its convergence without sacrificing accuracy.
.
mallela ankamma rao
2022-4-7
ok sir i will try sir.thanks alot for your support sir
with your help my project came last stage sir.
please help me in final work sir.with this my project will ending stage sir
how to plot graphs for simulation sir. i wrote this sir was it correct?
if not please modify sir
thank you very much for your help sir
hold on
figure(1),plot(t,x(:,1),'m--','linewidth',2) % for x(1)
xlabel('days');
ylabel('population');
hold on
figure(2),plot(t,x(:,2),'k--','linewidth',2) % for x(2)
xlabel('days');
ylabel('population');
hold on
figure(3),plot(t,x(:,3),'r--','linewidth',2) % for x(3)
xlabel('days');
ylabel('population');
hold on
figure(4),plot(t,x(:,4),'g--','linewidth',2) % for x(4)
xlabel('days');
ylabel('population');
hold on
figure(5),plot(t,x(:,5),'b--','linewidth',2) % for x(5)
xlabel('days');
ylabel('population');
hold on
figure(6),plot(t,x(:,6),'y--','linewidth',2) % for x(6)
xlabel('days');
ylabel('population');
hold on
figure(7),plot(t,x(:,7),'y--','linewidth',2) % for x(7)
xlabel('days');
ylabel('population');
hold on
figure(8),plot(t,x(:,8),'k','linewidth',2) % for x(8)
xlabel('days');
ylabel('population');
hold on
figure(9),plot(t,x(:,1),'m--','linewidth',2;t,x(:,2),'k--','linewidth',2;t,x(:,3),'r--','linewidth',2;t,x(:,4),'g--','linewidth',2;
t,x(:,5),'b--','linewidth',2;t,x(:,6),'y--','linewidth',2;t,x(:,7),'y--','linewidth',2;(t,x(:,8),'k','linewidth',2)
% for all x
plot (t,A_data,t,A) % for actual A_data and predicted A
A_data= [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
tspan = [1 2 3 4 5 7 9 11 13 15 17 21 25 30 40 50 60];
N = 1390000000;
e0 = 121737;
a0 = 150;
i0 = 100;
q0 = 30;
h0 = 10;
r0 = 10;
d0 = 0;
s0 = N-e0-a0-i0-q0-h0-r0-d0;
x0 = [s0; e0; a0; i0; q0; h0; r0; d0];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(A_data,tspan,x0),B0,lb,ub,options );
disp(B)
function A = immanuel(B,tspan,x0)
[T,Av] = ode45(@(t,x)DifEq(t, x),tspan, x0);
function dA = DifEq(t, x)
N = 1390000000;
pi = 70000;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,2);
end
mallela ankamma rao
2022-4-7
ok sir
.please see below code sir
please modify graph for simulation and graph for actual vs fit sir
thank you sir
A_data= [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
t = [1 2 3 4 5 7 9 11 13 15 17 21 25 30 40 50 60];
N = 1390000000;
e0 = 121737;
a0 = 150;
i0 = 100;
q0 = 30;
h0 = 10;
r0 = 10;
d0 = 0;
s0 = N-e0-a0-i0-q0-h0-r0-d0;
x0 = [s0; e0; a0; i0; q0; h0; r0; d0];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
B0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(B)immanuel(A_data,t,x0),B0,lb,ub,options );
disp(B)
function A = immanuel(B,t,x0)
[T,Av] = ode45(@(t,x)DifEq(t, x),t, x0);
plot (t,A_data,t,A) % for actual A_data and predicted A
hold on
figure(1),plot(time,Av(:,1),'m--','linewidth',2)
xlabel('days');
ylabel('population');
legend('s(t)');
title('susceptible')
hold on
figure(2),plot(time,Av(:,2),'k--','linewidth',2)
xlabel('days');
ylabel('population');
legend('I(t)');
title('insusceptible')
hold on
figure(3),plot(time,Cv(:,3),'r--','linewidth',2)
xlabel('days');
ylabel('population');
legend('A(t)');
title('Allied')
hold on
figure(4),plot(time,Av(:,4),'g--','linewidth',2)
xlabel('days');
ylabel('population');
legend('G(t)');
title('genered')
hold on
figure(5),plot(time,Av(:,5),'b--','linewidth',2)
xlabel('days');
ylabel('population');
legend('Q(t)');
title('Attempted')
hold on
figure(6),plot(time,Av(:,6),'y--','linewidth',2)
xlabel('days');
ylabel('population');
legend('U(t)');
title('unattempted')
hold on
figure(7),plot(time,Av(:,7),'y--','linewidth',2)
xlabel('days');
ylabel('population');
legend('R(t)');
title('Recovered')
hold on
figure(8),plot(time,Av(:,8),'k','linewidth',2)
xlabel('days');
ylabel('population');
legend('UN(t)');
title('unrecovered')
hold on
figure(9),plot(time,Av)
legend('Population')
xlabel('days')
ylabel('population')
title('maleria')
function dA = DifEq(t, x)
N = 1390000000;
pi = 70000;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = B(1);
alpha =B(2);
gamma = B(3);
upsilon =B(4);
epsilon = B(5);
lamda = B(6);
sigma = B(7);
kappa = B(8);
nu = B(9);
xi = B(10);
xdot = zeros(8,1);
xdot(1) = pi -beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zeta*x(3)+eta*x(4)+theta*x(5)+iota*x(6))*(x(1)/N) -(delta+mu)*x(2);
xdot(3) = rho*delta*x(2)-(lamda+gamma+nu+mu)*x(3);
xdot(4) = (1-rho)*delta*x(2)-(sigma+kappa+mu)*x(4);
xdot(5) = lamda*x(3) + sigma*x(4)-(alpha+upsilon+mu)*x(5);
xdot(6) = alpha*x(6) + kappa*x(4)- (epsilon+xi+mu)*x(6);
xdot(7) = gamma*x(3) + upsilon*x(5) + epsilon*x(6);
xdot(8) = nu*x(3) + xi*x(6);
dA = xdot;
end
A = Av(:,2);
end
Star Strider
2022-4-7
I did not run it again, however it appears to be correct.
The plot titles tell what the simulation is doing. More information is always better.
mallela ankamma rao
2022-4-7
figure(9),plot(time,Cv)
xlabel('days')
ylabel('population')
title('maleria')
i did not get all x(1),x(2),x(3),x(4,),x(5),x(6),x(7) and x(8) in this figure sir
for this simulation also running long time sir
figure,plot (time,c_data, time,Cv) % actual verses fitted
Also i did not get this figure sir
this two graphs only important to me sir
please modify this if necessary
the whole code is this sir if possible please run it
sir i requested you please spend 5 minutes for me and run the program and changes anything if need sir
this is correct code and final code i have sumbitted sir
thank you sir
code:
c_data= [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
time = [1 2 3 4 5 7 9 11 13 15 17 21 25 30 40 50 60];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
p0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) kinetics(p,time,c_data),p0,lb,ub,options );
disp(p)
function C=kinetics(p,time,c_data)
c0 = [1217378052,100,150,200,100,100,200,300];
[T,Cv]=ode45(@DifEq,time,c0);
hold on
figure,plot (time,c_data,time,Cv)
hold on
figure(1),plot(time,Cv(:,1),'m--','linewidth',2)
xlabel('days');
ylabel('population');
legend('s(t)');
title('susceptible')
hold on
figure(2),plot(time,Cv(:,2),'k--','linewidth',2)
xlabel('days');
ylabel('population');
legend('I(t)');
title('insusceptible')
hold on
figure(3),plot(time,Cv(:,3),'r--','linewidth',2)
xlabel('days');
ylabel('population');
legend('A(t)');
title('Allied')
hold on
figure(4),plot(time,Cv(:,4),'g--','linewidth',2)
xlabel('days');
ylabel('population');
legend('G(t)');
title('genered')
hold on
figure(5),plot(time,Cv(:,5),'b--','linewidth',2)
xlabel('days');
ylabel('population');
legend('Q(t)');
title('Attempted')
hold on
figure(6),plot(time,Cv(:,6),'y--','linewidth',2)
xlabel('days');
ylabel('population');
legend('U(t)');
title('unattempted')
hold on
figure(7),plot(time,Cv(:,7),'y--','linewidth',2)
xlabel('days');
ylabel('population');
legend('R(t)');
title('Recovered')
hold on
figure(8),plot(time,Cv(:,8),'k','linewidth',2)
xlabel('days');
ylabel('population');
legend('UN(t)');
title('unrecovered')
hold on
figure(9),plot(time,Cv)
legend('HEPITITES)
xlabel('days')
ylabel('population')
title('maleria')
function dC=DifEq(time,c)
N = 1390000000;
pi = 70000;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = p(1);
alpha =p(2);
gamma = p(3);
upsilon =p(4);
epsilon = p(5);
lamda = p(6);
sigma = p(7);
kappa = p(8);
nu = p(9);
xi = p(10);
dcdt = zeros(8,1);
dcdt(1) = pi -beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -mu*c(1);
dcdt(2) = beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -(delta+mu)*c(2);
dcdt(3) = rho*delta*c(2)-(lamda+gamma+nu+mu)*c(3);
dcdt(4) = (1-rho)*delta*c(2)-(sigma+kappa+mu)*c(4);
dcdt(5) = lamda*c(3) + sigma*c(4)-(alpha+upsilon+mu)*c(5);
dcdt(6) = alpha*c(6) + kappa*c(4)- (epsilon+xi+mu)*c(6);
dcdt(7) = gamma*c(3) + upsilon*c(5) + epsilon*c(6);
dcdt(8) = nu*c(3) + xi*c(6);
dC = dcdt;
end
C=c_data-Cv(:,2);
end
Torsten
2022-4-7
Why do you plot in "kinetics" and not in your script when lsqnonlin has finished ?
After you got the parameters from lsqnonlin, you can call "kinetics" to get Cv for your plot:
c_data= [503,502,575,513,484,376,405,406,339,490,444,397,341,397,341,356,387,359,346];
time = [1 2 3 4 5 7 9 11 13 15 17 21 25 30 40 50 60];
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
p0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) kinetics(p,time,c_data),p0,lb,ub,options );
c_data_minus_Cv=kinetics(p,time,c_data);
Cv = c_data - c_data_minus_Cv;
plot(time,c_data, time,Cv)
mallela ankamma rao
2022-4-13
good morning sir
i got parameter values from lsqnonlin. i got figure for c_data vs fitting .
i have doubt sir .was the graph is correct or not?
if not i request you please modify sir
thank you sir
my code;
time = data(:,1);
c_data= data(:,2);
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
p0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) immanuel(p,time,c_data),p0,lb,ub,options );
disp(p)
function C=immanuel(p,time,c_data)
c0 = [1217378052,100,10,5,3,3,1,1];
[T,Cv]=ode45(@DifEq,time,c0);
function dC=DifEq(time,c)
N = 1390000000;
pi = 150;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = p(1);
alpha =p(2);
gamma = p(3);
upsilon =p(4);
epsilon = p(5);
lamda = p(6);
sigma = p(7);
kappa = p(8);
nu = p(9);
xi = p(10);
dcdt = zeros(8,1);
dcdt(1) = pi -beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -mu*c(1);
dcdt(2) = beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -(delta+mu)*c(2);
dcdt(3) = rho*delta*c(2)-(lamda+gamma+nu+mu)*c(3);
dcdt(4) = (1-rho)*delta*c(2)-(sigma+kappa+mu)*c(4);
dcdt(5) = lamda*c(3) + sigma*c(4)-(alpha+upsilon+mu)*c(5);
dcdt(6) = alpha*c(6) + kappa*c(4)- (epsilon+xi+mu)*c(6);
dcdt(7) = gamma*c(3) + upsilon*c(5) + epsilon*c(6);
dcdt(8) = nu*c(3) + xi*c(6);
dC = dcdt;
end
C=c_data-Cv(:,2);
plot(time,c_data,time,C)
end
3 个评论
Torsten
2022-4-13
You mustn't plot
plot(time,c_data,time,C)
but
plot(time,c_data,time,Cv(:,2))
and - as I already wrote - you should do it in your calling program:
time = data(:,1);
c_data= data(:,2);
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
p0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) immanuel(p,time,c_data),p0,lb,ub,options );
Cdata_minus_Cv =immanuel(p,time,c_data);
Cv = c_data-Cdata_minus_Cv;
plot(time,[c_data,Cv])
mallela ankamma rao
2022-4-13
thanks alot for your reply dear Torsten sir
i tried sir but it showing error sir
i felt very difficult about graph sir
please resolve this sir
code:
data = [ 1 36456
2 70217
3 108526
4 145164
5 181275
6 213547
7 240654
8 272721
9 299078
10 333744
11 363778
12 394132
13 421468
14 443409
15 469810
16 487974
17 514736
18 541727
19 568561
20 593150
21 612324
22 636210
23 660446
24 684370
25 706720
26 725294
27 745627
28 761699
29 782228
30 804185
31 823231
32 842753
33 862213
34 881151
35 900118
36 918709
37 937090
38 955342
39 973350
40 991455
41 1009475
42 1026953
43 1044144
44 1061441
45 1077789
46 1093666
47 1109004
48 1124038
49 1138716
50 1153092
51 1167220
52 1181179
53 1195019
54 1208837
55 1222551
56 1236139
57 1249576
58 1262481
59 1276016
60 1289370
61 1302464
62 1315314
63 1328090
64 1340619
65 1353343
66 1365137
67 1376739
68 1388198
69 1399694
70 1411209
71 1422686
72 1434108
73 1445094
74 1456141
75 1467207
76 1478261
77 1489368
78 1500598
79 1511874
80 1523705
81 1535790
82 1548121
83 1560820
84 1573720
85 1586918
86 1600698
87 1614958
88 1629596
89 1644646
90 1659887
91 1675382
92 1691064
93 1706817
94 1722602
95 1738639
96 1754950
97 1771695
98 1788880
99 1806476
100 1824847];
time = data(:,1);
c_data= data(:,2);
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
p0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) immanuel(p,time,c_data),p0,lb,ub,options );
disp(p)
c_data_minus_Cv=immanuel(p,time,c_data);
Cv = c_data - c_data_minus_Cv;
plot(time,c_data, time,Cv)
c0 = [1217378052,100,10,5,3,3,1,1];
Cv = c_data - c_data_minus_Cv;
plot(time,c_data, time,Cv)
c0 = [1217378052,100,10,5,3,3,1,1];
[T,Cv]=ode45(@DifEq,time,c0);
function dC=DifEq(time,c)
N = 1390000000;
pi = 150;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = p(1);
alpha =p(2);
gamma = p(3);
upsilon =p(4);
epsilon = p(5);
lamda = p(6);
sigma = p(7);
kappa = p(8);
nu = p(9);
xi = p(10);
dcdt = zeros(8,1);
dcdt(1) = pi -beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -mu*c(1);
dcdt(2) = beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -(delta+mu)*c(2);
dcdt(3) = rho*delta*c(2)-(lamda+gamma+nu+mu)*c(3);
dcdt(4) = (1-rho)*delta*c(2)-(sigma+kappa+mu)*c(4);
dcdt(5) = lamda*c(3) + sigma*c(4)-(alpha+upsilon+mu)*c(5);
dcdt(6) = alpha*c(6) + kappa*c(4)- (epsilon+xi+mu)*c(6);
dcdt(7) = gamma*c(3) + upsilon*c(5) + epsilon*c(6)-mu*c(7);
dcdt(8) = nu*c(3) + xi*c(6);
dC = dcdt;
end
ERROR:
Unrecognized function or variable 'immanuel'.
Error in errrrr>@(p)kinetics(p,time,c_data) (line 20)
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) immanuel(p,time,c_data),p0,lb,ub,options );
Error in lsqnonlin (line 218)
initVals.F = feval(funfcn{3},xCurrent,varargin{:});
Error in errrrr (line 20)
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) immanuel(p,time,c_data),p0,lb,ub,options );
mallela ankamma rao
2022-4-13
dear torsten sir and dear star strider sir
i request you please reply anyone how to modify this code sir
how to overcome this error sir
Torsten
2022-4-13
data = [ 1 36456
2 70217
3 108526
4 145164
5 181275
6 213547
7 240654
8 272721
9 299078
10 333744
11 363778
12 394132
13 421468
14 443409
15 469810
16 487974
17 514736
18 541727
19 568561
20 593150
21 612324
22 636210
23 660446
24 684370
25 706720
26 725294
27 745627
28 761699
29 782228
30 804185
31 823231
32 842753
33 862213
34 881151
35 900118
36 918709
37 937090
38 955342
39 973350
40 991455
41 1009475
42 1026953
43 1044144
44 1061441
45 1077789
46 1093666
47 1109004
48 1124038
49 1138716
50 1153092
51 1167220
52 1181179
53 1195019
54 1208837
55 1222551
56 1236139
57 1249576
58 1262481
59 1276016
60 1289370
61 1302464
62 1315314
63 1328090
64 1340619
65 1353343
66 1365137
67 1376739
68 1388198
69 1399694
70 1411209
71 1422686
72 1434108
73 1445094
74 1456141
75 1467207
76 1478261
77 1489368
78 1500598
79 1511874
80 1523705
81 1535790
82 1548121
83 1560820
84 1573720
85 1586918
86 1600698
87 1614958
88 1629596
89 1644646
90 1659887
91 1675382
92 1691064
93 1706817
94 1722602
95 1738639
96 1754950
97 1771695
98 1788880
99 1806476
100 1824847];
time = data(:,1);
c_data= data(:,2);
beta0 = 0.5;
alpha0 = 0.004;
gamma0 = 0.1;
upsilon0 = 0.13;
epsilon0= 0.07;
lamda0= 0.1;
sigma0= 0.07;
kappa0= 0.03;
nu0 = 0.0001;
xi0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
p0 = [beta0; alpha0; gamma0; upsilon0; epsilon0; lamda0; sigma0; kappa0; nu0; xi0 ];
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@(p) immanuel(p,time,c_data),p0,lb,ub,options );
disp(p)
Cdata_minus_Cv =immanuel(p,time,c_data);
Cv = c_data-Cdata_minus_Cv;
plot(time,[c_data,Cv])
function C=immanuel(p,time,c_data)
c0 = [1217378052,100,10,5,3,3,1,1];
[T,Cv]=ode45(@DifEq,time,c0);
function dC=DifEq(time,c)
N = 1390000000;
pi = 150;
zeta = 0.1;
eta = 0.2;
theta = 0.3;
iota = 0.3;
delta = 0.1;
rho = 0.5;
mu = 0.0000425;
beta = p(1);
alpha =p(2);
gamma = p(3);
upsilon =p(4);
epsilon = p(5);
lamda = p(6);
sigma = p(7);
kappa = p(8);
nu = p(9);
xi = p(10);
dcdt = zeros(8,1);
dcdt(1) = pi -beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -mu*c(1);
dcdt(2) = beta*(zeta*c(3)+eta*c(4)+theta*c(5)+iota*c(6))*(c(1)/N) -(delta+mu)*c(2);
dcdt(3) = rho*delta*c(2)-(lamda+gamma+nu+mu)*c(3);
dcdt(4) = (1-rho)*delta*c(2)-(sigma+kappa+mu)*c(4);
dcdt(5) = lamda*c(3) + sigma*c(4)-(alpha+upsilon+mu)*c(5);
dcdt(6) = alpha*c(6) + kappa*c(4)- (epsilon+xi+mu)*c(6);
dcdt(7) = gamma*c(3) + upsilon*c(5) + epsilon*c(6);
dcdt(8) = nu*c(3) + xi*c(6);
dC = dcdt;
end
C=c_data-Cv(:,2);
end
14 个评论
mallela ankamma rao
2022-4-20
Good evening Torsten sir Good evening star strider sir How to draw contour plot two different set of values of beta={0.2,0.3} and gamma ={0.1,0.2}. Please give me any code relating this model sir
Thank you sir
Torsten
2022-4-20
How to draw contour plot two different set of values of beta={0.2,0.3} and gamma ={0.1,0.2}.
A contour plot of what ?
mallela ankamma rao
2022-4-20
A contour plot of reproduction number R = 1.5 with respect to rate beta = 0.2 and gamma = 0.1 Sir i have only this image sir But i dont know how to draw this please help me or give me some reference examples
mallela ankamma rao
2022-4-21
sir please tell me how to write this contour plot any suggestions or examples.
mallela ankamma rao
2022-4-22
good evening sir
sir i have to draw contour plot . i am taking some examples from other papers i am trying to find contour plot but i did not get results sir. please tell me where i had made mistake
the code is
% R0 = x*(1-b)*(1-v)*(u+y+p+t+w)*e/(u*(u+w)*(u+y+p+t))
% b =1.10340
% v =1.03330
% u =0.72936
% p =0.65722
% t =3.04270
% w =0.327688
% e =0.521749
% by substituting all values i got function Ro = x.*(0.00230).*((y+4.7569)/(y+4.442929))
x = linspace(0,0.02,1)
y = linspace(0,0.02,0.04)
[xm,ym] = meshgrid(x,y)
R0 = xm.*(0.00230).*((ym+4.7569)/(ym+4.442929))
contourf(xm,ym,R0)
colorbar
but original figure from other paper is
here alpha c is x and delta is y
i request you sir please help how to get this graph sir
thank you sir
Torsten
2022-4-22
R0 = xm.*0.00230.*(ym+4.7569)./(ym+4.442929)
instead of
R0 = xm.*(0.00230).*((ym+4.7569)/(ym+4.442929))
mallela ankamma rao
2022-4-22
sir one more doubt please tell me
the following code got graph for N=150 but i did not get graph for N=200 or 250
so please tell me sir how did i take N value.
code is
x = [0.125;0.2518;0.3518;0.4518;0.5518;0.6518;0.7518;0.8518;0.9518];
y = [0.0113;0.0213;0.0313;0.0413;0.0513;0.0613;0.0713;0.0813;0.0913];
z = [1.1033;1.6686;1.1801;1.6606;2.1217;2.5696;2.0085;2.4408;2.8680];
N = 250;
xv = linspace(min(x), max(x),N);
yv = linspace(min(y), max(y),N);
[X,Y] = meshgrid(xv, yv);
figure
contourf(X, Y, Z)
colorbar
xlabel('x')
ylabel('y')
error is
error using contourf (line 57)
The size of X must match the size of Z or the number of columns of Z.
Error in contour_plot (line 87)
contourf(X, Y, Z)
Torsten
2022-4-22
You must evaluate your function for Z at every element X(i,j),Y(i,j) such that X,Y and Z finally are of equal size.
mallela ankamma rao
2022-4-29
Good afternoon sir
sir please tell me how to draw graph for infected data vs predicted for this lscurvefit code.
i wrote code like this
plot(time,infected ,time,Cv) but it shows error sir.
code is
data = [ 1 25
2 75
3 227
4 296
5 258
6 236
7 192
8 126
9 71
10 28
11 11
12 7];
time = data(1:12,1);
infected= data(1:12,2);
beta0 = 1;
lamdaa0= 0.019;
lamdas0= 0.0715;
etas0 = 0.03;
etaq0 = 0.04;
gammaa0 = 0.2;
gammaq0 = 0.13;
gammah0 = 0.07;
mua0 = 0.0001;
muh0 = 0.0002;
lb =[0,0,0,0,0,0,0,0,0,0]; ub = [1,1,1,1,1,1,1,1,1,1];
B0 = [beta0; lamdaa0; lamdas0; etas0; etaq0; gammaa0; gammaq0; gammah0; mua0; muh0 ];
options=optimset('MaxFunEvals', 10000, 'MaxIter', 10000, 'TolFun', 0.0001, 'TolX',0.0001,'Display','on');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@diff1,B0,time,infected,lb,ub,options);
disp(B)
%plot(time,infected,time,C)
function C = diff1(B,time)
x0 = [1217378052,120,3,2,1,1,0,0];
[t,Cv] = ode45(@DifEq, time, x0);
function dC = DifEq(t, x)
N = 1390000000;
pi = 700 ;
zetaa = 0.1;
zetas = 0.2;
zetaq = 0.3;
zetah = 0.3;
omega = 0.2;
theta = 0.5;
mu = 0.0000425;
beta = B(1);
lamdaa= B(2);
lamdas = B(3);
etas = B(4);
etaq = B(5);
gammaa = B(6);
gammaq = B(7);
gammah = B(8);
mua = B(9);
muh = B(10);
xdot = zeros(7,1);
xdot(1) = pi -beta*(zetaa*x(3) +zetas*x(4) +zetaq*x(5)+zetah*x(6))*(x(1)/N) -mu*x(1);
xdot(2) = beta*(zetaa*x(3) +zetas*x(4) +zetaq*x(5)+zetah*x(6))*(x(1)/N) -(omega+mu)*x(2);
xdot(3) = theta*omega*x(2)-(lamdaa+gammaa+mua+mu)*x(3);
xdot(4) = (1-theta)*omega*E-(lamdas+etas+mu)*x(4);
xdot(5) = lamdaa*Ia+lamdas*Is-(etaq+gammaq+mu)*x(5);
xdot(6) = etas*Is+etaq*Q- (gammah+muh+mu)*x(6);
xdot(7) = gammaa*x(4) + gammaq*x(5) + gammah*x(6)
end
C = Cv(:,2);
plot(time,infected ,time,Cv)
end
Error is
Unrecognized function or variable 'infected'.
Error in lscurvfit6>diff1 (line 72)
plot(time,infected ,time,Cv)
Error in lsqcurvefit (line 225)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in lscurvfit6 (line 19)
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@diff1,B0,time,infected,lb,ub,options);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
mallela ankamma rao
2022-4-29
Mr torsten sir and Mr star strider sir please reply anyone
how can i draw graph for above code.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Acquisition Using Kinect for Windows Hardware 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)