estimation of parameters by using lsqnonlin for system of differential equations

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

[mallela ankamma rao]

sir please help anyone to solve this error

Answer 1

[Torsten]

[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.

[mallela ankamma rao]

thanks for your reply sir

please tell me how to modify this

[Torsten]

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]

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

Answer 2

[Star Strider]

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

.

[mallela ankamma rao]

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]

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]

@Star Strider

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]

thanks alot sir

i will try sir

[mallela ankamma rao]

[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]

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]

ok sir i will try . thank you somuch sir

[Star Strider]

My pleasure!

[mallela ankamma rao]

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]

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]

its worked sir.thanks alot Torston sir

and thanks alot star strider sir

[Star Strider]

Our pleasure!

[mallela ankamma rao]

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);

[Torsten]

A = Av(:,2);  % column of Av that corresponds to your measurement data

instead of

A = Av;

[Star Strider]

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]

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]

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]

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]

@mallela ankamma rao

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]

i got it . Thank you very much star strider sir

Thank you very much torsten sir.

[mallela ankamma rao]

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]

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]

ok sir i will try right now.

thanks alot sir.

[Star Strider]

My pleasure!

[mallela ankamma rao]

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]

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]

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]

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]

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]

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]

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

[Star Strider]

I would add a title to each plot tellling what the population is that is plotted.

[mallela ankamma rao]

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]

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]

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]

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]

ok sir i can try

thanks alot sir

Answer 3

[mallela ankamma rao]

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

[Torsten]

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]

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]

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

Answer 4

[Torsten]

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

[mallela ankamma rao]

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]

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]

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]

sir please tell me how to write this contour plot any suggestions or examples.

[mallela ankamma rao]

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]

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]

thanks alot for your reply torsten sir.

[mallela ankamma rao]

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]

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]

Thanks alot sir I requesting you please provide me any examples sir

[Torsten]

Do you have a formula how the z-vector is computed as a function of x and y ?

[mallela ankamma rao]

Thanks alot sir I am preparing the formula sir.

[mallela ankamma rao]

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]

Mr torsten sir and Mr star strider sir please reply anyone

how can i draw graph for above code.

[Pavl M.]

Dear Friend,

How are you?

Are you still on this, what is the approximated volume of research works requests/clients/order?

There are already reliable, very valued original few real world mathematical and physical models and ethical framework lawful business proposal in our disposal. Let's make from it "limonade" furthermore.

Let's make order in it. Hire me continuously, stable 1 here or around for good order in it.

I havn't seen your works from before, so it is new to me and I just do it to correct, amend, order, fix:

Let's understand firstly what it factually is, the underlying mechanics, principles of correct operation and right directed evolution of the system and processes towards product only.

So far which biological or chemical process it is all about?

Where is original

E =;

Ia = ;

Is = ;

Q = ;

?

Who I am? Please get to know better, all datum about my honest work and utilization are provided via http links posted here and there.

We can work hybrid as well, outdoor and indoor (demonstrated already kept connected, kindly don't dissapoint, hire , wait, get).

I need investement to upgrade to enterprise service hardware and also to reward so valuable time, not into software and text complications. Kindly command/order them to invest to me, not others disturbing, wha you, invest you if you will be happy that I will continue to serve your objectives normally well.

Let's get on, Pitch Deck stakeholders, make necessary remittance transactions for good job we really did before and go ahead continue doing this.

clc
clear all
close all

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', 1100, 'MaxIter', 1100, 'TolFun', 0.00001, 'TolX',0.00001,'Display','on');
[B,resnorm,RESIDUAL,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqcurvefit(@diff1,B0,time,infected,lb,ub,options);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

disp(B)

    1.0000
    0.0000
    0.0000
    0.0000
    1.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000

OUTPUT

OUTPUT = struct with fields:
      firstorderopt: 4.6276e-07
         iterations: 5
          funcCount: 66
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 0.0032
            message: 'Local minimum found....'
       bestfeasible: []
    constrviolation: []

C = RESIDUAL

C = 12×1

95.0000 24.7749 -142.8381 -223.7400 -194.6765 -179.2484 -139.9337 -77.1106 -24.0776 17.9305

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(time,infected,time,C)

function C = diff1(B,time)

x0 =  [1217378052,120,3,2,1,1,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);
        %ToDo: update the parameters
        E =0.521749;
        Ia = 1;
        Is = 1;
        Q = 1;
        
        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);
        dC = xdot;
      
  end
 C = Cv(:,2);
end