Unable to perform assignment because the left and right sides have a different number of elements.

3 次查看（过去 30 天）
Fares
Fares 2024-2-9
编辑： Fares 2024-2-9
Ran in:
I am trying to run this code but always get this error
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;
eta=0.05;
sigma=5;
zeta=0.05;
omega0=0.001;
NN1=5000;
NN2=5000;
TT=linspace(0.01,350,NN1);
syms M b Msol positive
Nut(M) = mu+(rho*M/(1+M)); % N(M)
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M))); % r(M)
lam(M) = 1/(1+Nut(M)); % λ(Μ)
P(M) = epsilon*M/(1+M);
dNut(M) = diff(Nut(M),M); % Ν'(Μ)
dgro(M) = diff(gro(M),M); % r'(M)
dlam(M) = diff(lam(M),M); % λ'(M)
dP(M) = diff(P(M),M); % P'(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));
b_num = linspace(0.001,8,NN2);
for u = 1:numel(b_num)
Msol(u) = vpa(solve(subs(num,b,b_num(u))));
end
Unable to perform assignment because the left and right sides have a different number of elements.

Error in sym/privsubsasgn (line 1200)
L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in indexing (line 1031)
C = privsubsasgn(L,R,inds{:});
However when the vector b_num starts from 0.008, the code works. Is there a way I can include the values of b_num less than 0.008? Many thanks!

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-2-9
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;
eta=0.05;
sigma=5;
zeta=0.05;
omega0=0.001;
NN1=5000;
NN2=5000;
TT=linspace(0.01,350,NN1);
syms M b Msol positive
Nut(M) = mu+(rho*M/(1+M)); % N(M)
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M))); % r(M)
lam(M) = 1/(1+Nut(M)); % λ(Μ)
P(M) = epsilon*M/(1+M);
dNut(M) = diff(Nut(M),M); % Ν'(Μ)
dgro(M) = diff(gro(M),M); % r'(M)
dlam(M) = diff(lam(M),M); % λ'(M)
dP(M) = diff(P(M),M); % P'(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));
b_num = linspace(0.001,8,NN2);
Msol = nan(size(b_num));
for u = 1:numel(b_num)
sol = vpa(solve(subs(num,b,b_num(u))));
if ~isempty(sol)
Msol(u) = sol;
end
end
  3 个评论
显示 1更早的评论
Torsten
Torsten 2024-2-9
编辑：Torsten 2024-2-9
Ran in:
What is the old one ? I thought it did not work.
If you want a fast solution, you must use the numerical function "roots" to solve for the roots of your polynomial in M of degree 7 instead of the symbolic "solve".
If your equation has more than one positive solution for M given b, I took the first one.
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;
eta=0.05;
sigma=5;
zeta=0.05;
omega0=0.001;
NN1=5000;
NN2=5000;
TT=linspace(0.01,350,NN1);
syms M b
Nut(M) = mu+(rho*M/(1+M)); % N(M)
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M))); % r(M)
lam(M) = 1/(1+Nut(M)); % λ(Μ)
P(M) = epsilon*M/(1+M);
dNut(M) = diff(Nut(M),M); % Ν'(Μ)
dgro(M) = diff(gro(M),M); % r'(M)
dlam(M) = diff(lam(M),M); % λ'(M)
dP(M) = diff(P(M),M); % P'(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));
p = fliplr(coeffs(num,M));
p = matlabFunction(p);
b_num = linspace(0.001,8,NN2);
Msol = nan(size(b_num));
for u = 1:numel(b_num)
sol = roots(p(b_num(u)));
sol = sol(real(sol)>0 & abs(imag(sol)) < 1e-8);
if ~isempty(sol)
Msol(u) = sol(1);
end
end
plot(b_num,Msol)
grid on

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2024-2-9
编辑：Matt J 2024-2-9
Ran in:
Your code assumes that solve() will always return one and only one solution. There is no reason to think that will always be the case (see below). So, the question becomes, what do you want to do when situations like in the examples below arise?
syms x real
vpa(solve(x^2==1))
ans = 
vpa(solve(x^2==-1))
ans = Empty sym: 0-by-1
  5 个评论
显示 3更早的评论
Torsten
Torsten 2024-2-9
编辑：Torsten 2024-2-9
By removing the semicolon behind the line
sol = roots(p(b_num(u)));
you can see the numerical values of the 7 computed roots before the negative and complex ones are rejected by the line
sol = sol(real(sol)>0 & abs(imag(sol)) < 1e-8);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by