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

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!

Answer 1

[Torsten]

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

[Fares]

Thank you Torsten. One more thing: this code takes longer than the old one. Is there a way to make it faster?

[Torsten]

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

[Fares]

Thanks a lot Torsten! you've been very helpful.

Answer 2

[Matt J]

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
 

[Fares]

Yes, Matt. I want the code to give me only one positive root of the equation for each value of b. You mean for b<0.008, there might be no positive root or there might be a complex root?

[Matt J]

Or multiple roots... We don't know what criteria you would apply to reject a susbet of the solutions in that case.

[Fares]

Thank you Matt for the clarification. That is very helpful!

[Torsten]

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

[Fares]

That is a good idea to keep in mind. Thanks Torsten for your continuous support!