Roots of an equation when two parameters are changed

1 次查看(过去 30 天)
I have this equation
r0=0.05;
k1=0.5;
k2=0.5;
mu=0.5;
rho=0.5;
epsilon=0.25;
K=1;
alpha=0.1;
q=0.1;
b=0.8;
zeta=0.075;
omega0=0.001
syms M sigma eta Msol positive
Nut(M) = mu+(rho*M/(1+M));
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M)));
lam(M) = 1/(1+Nut(M));
P(M) = epsilon*M/(1+M);
eqn = (q./b).*gro(M).*(eta+P(M)).*(1+Nut(M)).*(1-(((alpha.*sigma.*q.*(eta+P(M)).*(1+Nut(M))+b.*(q+alpha).*(zeta.*M-omega0)))./(b.*alpha.*sigma.*K)))-(q./sigma).*(zeta.*M-omega0)==0;
I can compute the roots of this equation when sigma is varied and eta = 0.05 using the command
[num,den]=numden(lhs(eqn));
sigma_num = linspace(0.01,1,20);
for u = 1:numel(sigma_num)
Msol(u) = vpa(solve(subs(num,sigma,sigma_num(u))));
end
which will produce 20 roots, each one is related to each value of sigma. Now, I would like to compute the roots of the equation when two parameters are varied, say sigma and eta such that eta_num = linspace(0.01,1,20) as well. I believe I should end up having 400 roots but How I can do that? Any help is appreciated! Many thanks!

采纳的回答

Torsten
Torsten 2023-6-19
编辑:Torsten 2023-6-19
"num" is a polynomial of degree 7 in M. Thus it has seven roots. But your code only gives one of these seven. So I hope you get the single root out of seven that you are after.
r0=0.05;
k1=0.5;
k2=0.5;
mu=0.5;
rho=0.5;
epsilon=0.25;
K=1;
alpha=0.1;
q=0.1;
b=0.8;
zeta=0.075;
omega0=0.001;
syms M sigma eta Msol positive
Nut(M) = mu+(rho*M/(1+M));
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M)));
lam(M) = 1/(1+Nut(M));
P(M) = epsilon*M/(1+M);
eqn = (q./b).*gro(M).*(eta+P(M)).*(1+Nut(M)).*(1-(((alpha.*sigma.*q.*(eta+P(M)).*(1+Nut(M))+b.*(q+alpha).*(zeta.*M-omega0)))./(b.*alpha.*sigma.*K)))-(q./sigma).*(zeta.*M-omega0)==0;
[num,den]=numden(lhs(eqn));
sigma_num = linspace(0.01,1,20);
eta_num = linspace(0.01,1,20);
n = numel(eta_num);
m = numel(sigma_num);
M_num = zeros(n,m);
%M_num = cell(n,m);
for i = 1:n
for j = 1:m
M_num(i,j) = vpa(solve(subs(num,[eta sigma],[eta_num(i),sigma_num(j)])==0,M));
%M_num{i,j} = roots(sym2poly(subs(num,[eta sigma],[eta_num(i),sigma_num(j)])));
end
end
M_num
M_num = 20×20
0.0135 0.0146 0.0157 0.0168 0.0180 0.0193 0.0206 0.0219 0.0233 0.0248 0.0263 0.0279 0.0296 0.0314 0.0333 0.0352 0.0373 0.0395 0.0418 0.0442 0.0143 0.0194 0.0247 0.0302 0.0360 0.0420 0.0482 0.0547 0.0615 0.0686 0.0760 0.0837 0.0918 0.1002 0.1090 0.1181 0.1276 0.1376 0.1479 0.1587 0.0150 0.0240 0.0333 0.0431 0.0532 0.0638 0.0748 0.0862 0.0982 0.1106 0.1235 0.1369 0.1509 0.1654 0.1804 0.1959 0.2120 0.2286 0.2458 0.2635 0.0157 0.0284 0.0416 0.0554 0.0697 0.0847 0.1003 0.1164 0.1332 0.1507 0.1687 0.1874 0.2068 0.2268 0.2474 0.2687 0.2905 0.3130 0.3360 0.3597 0.0164 0.0326 0.0495 0.0672 0.0856 0.1047 0.1246 0.1452 0.1666 0.1887 0.2116 0.2352 0.2595 0.2845 0.3102 0.3366 0.3636 0.3912 0.4194 0.4481 0.0170 0.0366 0.0571 0.0785 0.1007 0.1238 0.1478 0.1727 0.1983 0.2249 0.2522 0.2803 0.3091 0.3387 0.3690 0.4000 0.4315 0.4637 0.4965 0.5297 0.0176 0.0405 0.0644 0.0893 0.1152 0.1421 0.1699 0.1987 0.2284 0.2590 0.2905 0.3227 0.3558 0.3895 0.4240 0.4591 0.4948 0.5311 0.5679 0.6053 0.0182 0.0442 0.0713 0.0996 0.1289 0.1594 0.1909 0.2234 0.2569 0.2913 0.3266 0.3627 0.3995 0.4371 0.4754 0.5142 0.5537 0.5937 0.6342 0.6752 0.0187 0.0477 0.0779 0.1094 0.1420 0.1759 0.2108 0.2468 0.2838 0.3217 0.3605 0.4001 0.4405 0.4816 0.5233 0.5656 0.6085 0.6519 0.6958 0.7401 0.0193 0.0510 0.0842 0.1187 0.1545 0.1915 0.2296 0.2689 0.3091 0.3503 0.3924 0.4353 0.4789 0.5232 0.5681 0.6135 0.6595 0.7060 0.7529 0.8002

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Mathematics 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by