Stability analysis of a non-linear ODE system
25 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I solved the following ODE system using the code:
syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta theta alpha CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);
[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]=solve(e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P)
with stability points for C_e = mu/theta, C_p = (k6*mu - alpha*theta + k12*mu)/(k5*theta), C_f = -(alpha*theta - k12*mu)/(k4*theta) and E = (k7*mu)/(k8*theta). All other variables are dependent on an input variable labelled CL. I want to do a stability analysis around these points.
I've determined the Jacobian using the code:
syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta theta alpha CL
f = [mu - eta*Sci*C,theta*Ce - eta*Sci*C,k1*Sci - k2*Sr,k1*Sci - k10*Sh,p1*Sr - k11*R - k14*P*R,k3*CL*R - k4*Cf,k4*Cf - k5*Cp + k6*Ce,k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E,k7*Ce - k8*E,p2*Sh - k15*HR*Ce,alpha - k9*HR*H,k1*Sci - k13*Sp,p3*Sp - k16*P];
J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
J = simplify(J)
I now have to determine the eigenvalues by first substituting equilibrium point values for the variables in the Jacobian with the steady state values ( in terms of CL) generated in the code above. I'm not sure how to do this in Matlab?
0 commentaires
Réponse acceptée
Torsten
le 7 Avr 2023
Modifié(e) : Torsten
le 7 Avr 2023
You won't succeed in this generality. To determine the eigenvalues, MATLAB had to solve for the roots of a polynomial of degree 13 with symbolic coefficients. This is in general only possible for polynomials up to degree 4. So you have to give values to the parameters of your function, I guess.
syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta theta alpha CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);
[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq]=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13],[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
f = [rhs(e1),rhs(e2),rhs(e3),rhs(e4),rhs(e5),rhs(e6),rhs(e7),rhs(e8),rhs(e9),rhs(e10),rhs(e11),rhs(e12),rhs(e13)];
%f = [mu - eta*Sci*C,theta*Ce - eta*Sci*C,k1*Sci - k2*Sr,k1*Sci - k10*Sh,p1*Sr - k11*R - k14*P*R,k3*CL*R - k4*Cf,k4*Cf - k5*Cp + k6*Ce,k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E,k7*Ce - k8*E,p2*Sh - k15*HR*Ce,alpha - k9*HR*H,k1*Sci - k13*Sp,p3*Sp - k16*P];
J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
J = simplify(J)
Jeq = subs(J,[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P],[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq])
eig(Jeq)
2 commentaires
Plus de réponses (1)
Torsten
le 13 Avr 2023
Modifié(e) : Torsten
le 13 Avr 2023
syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta alpha theta CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);
%Now solve to find the steady state values in terms of (A, B, C, D, Cl and the parameters as listed), then determine the Jacobian J:
[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq]=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13],[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
f = [rhs(e1),rhs(e2),rhs(e3),rhs(e4),rhs(e5),rhs(e6),rhs(e7),rhs(e8),rhs(e9),rhs(e10),rhs(e11),rhs(e12),rhs(e13)];
J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
J = simplify(J)
%Substitute the steady state values obtained into the Jacobian:
Jeq = subs(J,[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P],[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq])
%Lower (lb) and upper (ub) bounds for the free parameters
lbk1 = 3-0.5*3;
lbk2 = 2-0.5*2;
lbk3 = 5-0.5*5;
lbk4 = 4-0.5*4;
lbk5 = 5-0.5*5;
lbk6 = 1-0.5*1;
lbk7 = 4-0.5*4;
lbk8 = 3-0.5*3;
lbk9 = 1-0.5*1;
lbk10 = 1.2-0.5*1.2;
lbk11 = 1-0.5*1;
lbk12 = 9-0.5*9;
lbk13 = 1-0.5*1;
lbk14 = 1-0.5*1;
lbk15 = 1-0.5*1;
lbk16 = 1-0.5*1;
lbp1 = 30-0.5*30;
lbp2 = 4-0.5*4;
lbp3 = 1.5-0.5*1.5;
lbmu = 25-0.5*25;
lbeta = 10-0.5*10;
lbalpha = 10-0.5*10;
lbtheta = 10-0.5*10;
lbCL = 0.5-0.5*0.5;
ubk1 = 3+0.5*3;
ubk2 = 2+0.5*2;
ubk3 = 5+0.5*5;
ubk4 = 4+0.5*4;
ubk5 = 5+0.5*5;
ubk6 = 1+0.5*1;
ubk7 = 4+0.5*4;
ubk8 = 3+0.5*3;
ubk9 = 1+0.5*1;
ubk10 = 1.2+0.5*1.2;
ubk11 = 1+0.5*1;
ubk12 = 9+0.5*9;
ubk13 = 1+0.5*1;
ubk14 = 1+0.5*1;
ubk15 = 1+0.5*1;
ubk16 = 1+0.5*1;
ubp1 = 30+0.5*30;
ubp2 = 4+0.5*4;
ubp3 = 1.5+0.5*1.5;
ubmu = 25+0.5*25;
ubeta = 10+0.5*10;
ubalpha = 10+0.5*10;
ubtheta = 10+0.5*10;
ubCL = 0.5+0.5*0.5;
%Number of trials
n = 100;
found = zeros(n,1);
% Make a Monte Carlo simulation
for i = 1:n
% Random values for the parameters uniformly distributed between their
% lower and upper bounds
k1num(i) = lbk1 + rand()*(ubk1-lbk1);
k2num(i) = lbk2 + rand()*(ubk2-lbk2);
k3num(i) = lbk3 + rand()*(ubk3-lbk3);
k4num(i) = lbk4 + rand()*(ubk4-lbk4);
k5num(i) = lbk5 + rand()*(ubk5-lbk5);
k6num(i) = lbk6 + rand()*(ubk6-lbk6);
k7num(i) = lbk7 + rand()*(ubk7-lbk7);
k8num(i) = lbk8 + rand()*(ubk8-lbk8);
k9num(i) = lbk9 + rand()*(ubk9-lbk9);
k10num(i) = lbk10 + rand()*(ubk10-lbk10);
k11num(i) = lbk11 + rand()*(ubk11-lbk11);
k12num(i) = lbk12 + rand()*(ubk12-lbk12);
k13num(i) = lbk13 + rand()*(ubk13-lbk13);
k14num(i) = lbk14 + rand()*(ubk14-lbk14);
k15num(i) = lbk15 + rand()*(ubk15-lbk15);
k16num(i) = lbk16 + rand()*(ubk16-lbk16);
p1num(i) = lbp1 + rand()*(ubp1-lbp1);
p2num(i) = lbp2 + rand()*(ubp2-lbp2);
p3num(i) = lbp3 + rand()*(ubp3-lbp3);
munum(i) = lbmu + rand()*(ubmu-lbmu);
etanum(i) = lbeta + rand()*(ubeta-lbeta);
alphanum(i) = lbalpha + rand()*(ubalpha-lbalpha);
thetanum(i) = lbtheta + rand()*(ubtheta-lbtheta);
CLnum(i) = lbCL + rand()*(ubCL-lbCL);
Jeqnum = subs(Jeq,[k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta alpha theta CL],[k1num(i) k2num(i) k3num(i) k4num(i) k5num(i) k6num(i) k7num(i) k8num(i) k9num(i) k10num(i) k11num(i) k12num(i) k13num(i) k14num(i) k15num(i) k16num(i) p1num(i) p2num(i) p3num(i) munum(i) etanum(i) alphanum(i) thetanum(i) CLnum(i)]);
% Compute eigenvalues
v{i} = double(eig(Jeqnum));
% Check if all real parts are <= 0
if any(real(v{i})>0)
found(i) = 1;
end
end
% List the number of trials where eigenvalues with real part > 0 occured
sum(found)
8 commentaires
Torsten
le 23 Avr 2023
Modifié(e) : Torsten
le 23 Avr 2023
Your output of parameters is after the loop over i, thus for i = n.
If found(n) = 1 by chance, you got a set of parameters for which at least one eigenvalue of the Jacobian is > 0.
But this will usually not be the case.
Try to understand the code and why you have to output the parameters here:
if any(real(v{i})>0)
found(i) = 1;
k1num(i),k2num(i),....
end
Voir également
Catégories
En savoir plus sur Equation Solving dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!