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

1 view (last 30 days)
Hi.I have a code error.
error:Unable to perform assignment because the left and right sides have a different number of elements.
Error in rate_eq_program_1 (line 26)
ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);
%%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%
q = 1.6021766208e-19;
c = 2.99792458e8;
kk = 1.38e-23;
kk_eV = kk/q;
T = 300;
h = 6.626068e-34;
h_eV = h/q;
hbar = h/(2*pi);
hbar_eV = hbar/q;
epsilon0 = 8.854187817e-12;
P = -1;
rL = 1;
rB = 0.2 ;
tB = sqrt(1-rB^2);
Qt =500;
% Qc =180;
Qi = 14300;
Qp = 10000;
%Qc = 1800;
L = 5e-6;
ng = 3.5;
sigma = 3;
m = 25;
A = 0.105e-12;
confine = 0.5;
Vp = 0.576e-18;
Vc = confine*A*L;
gN =5e-20;
Ntr = 1e24;
tc = 0.5e-9;
%tc = 0.7e-9;
alphai = 1000;
%alpha = 1.5;
alpha =1;
confinec = 0.3;
vg = c/ng;
vgp = c/6/ng;
Vpc = 0.22e-18;
KD = confinec*vgp*alpha*gN/2;
KA = -confinec*vgp*gN/2;
lambda0 = 1554e-9;
deltalambda = 0.01e-9;
lambda_start = 1540e-9;
lambda_end = 1560e-9;
lambda = lambda_start:deltalambda:lambda_end;
omega = 2*pi*c./lambda;
omega0 = omega(find(lambda==lambda0));
E = hbar_eV.*omega;
gamai = omega./2/Qi;
gamap = omega./2/Qp;
gamat = omega./2/Qt;
%gamac = omega./2/Qc;
%gamac = omega./2/Qc;
gamac = gamat-gamai-gamap;
gama1 = gamac./2;
gama2 = gama1;
COS2 = 0.5*gama1./gama1.*tB^2/rB-0.5*tB^2/rB-rB;
SIN2 = P*tB./(2.*gama1*rB).*sqrt(4.*gama1.*gama1 - tB^2.*(gamac.^2));
% CS12 = 1/(1i*tB).*(1i*SIN2 + rB);
GLA = gN*vgp./(E.*Vpc*q);
GN = confine*vg*gN;
tin = 2*L/vg;
M = 0.5*(length(E)-1);
Tmin = 0;
Tmax = 5e-9;
dt = 0.01e-12;
Time = Tmin:dt:Tmax;
N0 = zeros(1,4);
N0(3) = 0.2;
%I =2e-3;TA I =12e-3;
I =0.5e-3;
Vnc= 0.24e-18;
%%%%%%%%%%%%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%%%%%%%%%%%
Time = Tmin:dt:Tmax;
z0 = [0.2,0.2,0.2,0.2];
[t,z] = ode45(@rate_eq_program_1,Time,z0);
param_rate_eq
figure(1)
plot(t,z(:,1)); % divided to normalize
xlabel('time [ns]','FontSize',14); % size of x label
ylabel('Arbitrary units','FontSize',14); % size of y label
set(gca,'FontSize',14); % size of tick marks on both axis
legend('N', 'Location','SE') % legend inside the plot
figure(2)
plot(t,z(:,2)); % divided to normalize
xlabel('time [ns]','FontSize',14); % size of x label
ylabel('Arbitrary units','FontSize',14); % size of y label
set(gca,'FontSize',14); % size of tick marks on both axis
legend('NC', 'Location','SE') % legend inside the plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function %%%%%%%%%%%%%%%%%%%%%%%%
function ys = rate_eq_program_1(t,z)
param_rate_eq_fano
ys=zeros(4,1);
%rR=(-p*gamma_c)/(1i*(delta_c)+gamma_T-1/2*(1-1i*henry)*conf_NC*V_g*g_n*(z(2)-N_0));
deltawc=KD.*(z(2))-1i*KA.*(z(2));
deltac=omega0+deltawc-omega;
rR=(-P.*gamac)./(1i*(deltac)+gamat-(1/2*(1-1i)*confinec*vg*gN).*(z(2)-Ntr));
rRS=-P.*gamac./(1i.*(deltac)+gamat);
%rRS=rB+(-1i*tB-rB).*gamac./(1i.*(deltac)+gamat);
g=alphai+(1/(2*L)).*log(1./(abs(rRS)).^2);
A=(2*epsilon0*n*ng)./(hbar.*omega0);
B=(1+(abs(rRS))).*(1-(abs(rRS)));
D=((g)-alphai);
H=c./(omega0.*n);
F=imag(rRS)./abs(rRS);
sigmas=A.*((B./D)+(H.*F));
ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);
ys(2) = (-z(2)./tc)-GLA.*(z(2)-Ntr).*(abs(z(4)).^2)./Vnc;
ys(3) = 1/2*(1-1i)*GN.*(z(1)-Nss).*z(3)+(1/tin.*((z(4)./rR)-z(3)));
ys(4) =((-1i.*deltawc-gamat).*z(4))-(P.*gamac.*z(3))+(1/2*(1-1i).*confinec*vg*gN.*(z(2)-Ntr).*z(4));
ys = ys';
end
i tried ys=zeros(4,M) But it's also a mistake.
i tried ys=zeros(4,length(E)) But it's also a mistake.
....
Thanks for helping someone.
  6 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 1 Feb 2020
Your function rate_eq_program_1 uses a number of variables defined in the main program, but without importing them. You need to add a "function" statement at the very top of your code, and add an "end" statement at the end of your code, in order to make rate_eq_program_1 into a nested function that can access the variables.
Now, in the main part of your program you have
lambda = lambda_start:deltalambda:lambda_end;
so lambda is a vector
omega = 2*pi*c./lambda;
so omega is a vector because it is built from the vector lambda
omega0 = omega(find(lambda==lambda0));
As mentioned before due to floating point round-off the == might not find any matches. It will find either 0 or 1 match, so omega0 will end up either empty or a scalar.
E = hbar_eV.*omega;
E is constructed from the vector omega, so E is a vector. Did you want to construct E from omega0 ?
Now, inside rate_eq_program_1 you have
ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);
As noted above, E is a vector, so the right hand side of the assignment is a vector, but the left hand side only has room for a scalar.
Now look a little further and it turns out that sigmas is a vector as well:
deltac=omega0+deltawc-omega;
You use all of the vector omega, so deltac will be a vector.
rR=(-P.*gamac)./(1i*(deltac)+gamat-(1/2*(1-1i)*confinec*vg*gN).*(z(2)-Ntr));
rRS=-P.*gamac./(1i.*(deltac)+gamat);
Those use the vector deltac so they will both be vectors.
g=alphai+(1/(2*L)).*log(1./(abs(rRS)).^2);
Uses the vector rRS so g is a vector
A=(2*epsilon0*n*ng)./(hbar.*omega0);
Use the undefined variable n . If we assume that n is a scalar, then A will be a scalar.
B=(1+(abs(rRS))).*(1-(abs(rRS)));
Uses the vector rRS so B will be a vector
D=((g)-alphai);
g is a vector so D will be a vector.
H=c./(omega0.*n);
Use the undefined variable n. If we assume that n is a scalar, then H will be a scalar.
F=imag(rRS)./abs(rRS);
Uses the vector rRS so F will be a vector.
sigmas=A.*((B./D)+(H.*F));
Uses the vector B, the vector D, the vector F, so sigmas will be a vector.
ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);
Uses the vector sigmas, so even if E was a scalar instead of being a vector, the right hand side will be a vector.
Now, if E were a scalar instead of a vector, and if you could make deltac a scalar instead of a vector, then the other variables involved would become scalars and the assignment could work.
ys(2) = (-z(2)./tc)-GLA.*(z(2)-Ntr).*(abs(z(4)).^2)./Vnc;
GLA is a vector constructed in the main routine based on E, so if E were a scalar then GLA would be a scalar.
ys(3) = 1/2*(1-1i)*GN.*(z(1)-Nss).*z(3)+(1/tin.*((z(4)./rR)-z(3)));
Uses the undefined variable Nss. Uses rR which as explained above is a vector because it is built from the vector deltac so if deltac were made into a scalar somehow, then rR would become a scalar and the right hand side would be a scalar (assuming that the undefined Nss is a scalar.)
ys(4) =((-1i.*deltawc-gamat).*z(4))-(P.*gamac.*z(3))+(1/2*(1-1i).*confinec*vg*gN.*(z(2)-Ntr).*z(4));
This used the vectors gamat and gammac from the main routine.
gamai = omega./2/Qi;
gamap = omega./2/Qp;
gamat = omega./2/Qt;
Those use all of the vector omega, so all three are vectors.
gamac = gamat-gamai-gamap;
The right hand side are all vectors, so gamac is a vector.
So ys(4) is a vector.
Beyond defining the undefined n, Nss, and param_rate_eq_fano, there are two approaches you could use:
  1. Loop the whole ode45 call over individual values implied by omega, such as by passing in an index and indexing with that were appropriate so that you get 4 x 1 output from the function; Or
  2. Reconstruct the ode45 not as a 4 x 1 system but instead as a (4*2001) x 1 system, running all of those 2001 odes at the same time.
  3 Comments

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 2 Feb 2020
The plotting needs to be fixed on this. And it takes a long time!!
function rate_eq_program
%%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%
q = 1.6021766208e-19;
c = 2.99792458e8;
kk = 1.38e-23;
kk_eV = kk/q;
T = 300;
h = 6.626068e-34;
h_eV = h/q;
hbar = h/(2*pi);
hbar_eV = hbar/q;
epsilon0 = 8.854187817e-12;
P = -1;
rL = 1;
rB = 0.2 ;
tB = sqrt(1-rB^2);
Qt =500;
% Qc =180;
Qi = 14300;
Qp = 10000;
%Qc = 1800;
L = 5e-6;
ng = 3.5;
sigma = 3;
m = 25;
A = 0.105e-12;
confine = 0.5;
Vp = 0.576e-18;
Vc = confine*A*L;
gN =5e-20;
Ntr = 1e24;
tc = 0.5e-9;
%tc = 0.7e-9;
alphai = 1000;
%alpha = 1.5;
alpha =1;
confinec = 0.3;
vg = c/ng;
vgp = c/6/ng;
Vpc = 0.22e-18;
KD = confinec*vgp*alpha*gN/2;
KA = -confinec*vgp*gN/2;
lambda0 = 1554e-9;
deltalambda = 0.01e-9;
lambda_start = 1540e-9;
lambda_end = 1560e-9;
lambda = lambda_start:deltalambda:lambda_end;
omega = 2*pi*c./lambda;
omega0 = omega(lambda==lambda0); %MARK
E = hbar_eV.*omega;
gamai = omega./2/Qi;
gamap = omega./2/Qp;
gamat = omega./2/Qt;
%gamac = omega./2/Qc;
%gamac = omega./2/Qc;
gamac = gamat-gamai-gamap;
gama1 = gamac./2;
gama2 = gama1;
COS2 = 0.5*gama1./gama1.*tB^2/rB-0.5*tB^2/rB-rB;
SIN2 = P*tB./(2.*gama1*rB).*sqrt(4.*gama1.*gama1 - tB^2.*(gamac.^2));
% CS12 = 1/(1i*tB).*(1i*SIN2 + rB);
GLA = gN*vgp./(E.*Vpc*q);
GN = confine*vg*gN;
tin = 2*L/vg;
M = 0.5*(length(E)-1);
Tmin = 0;
Tmax = 5e-9;
dt = 0.01e-12;
Time = Tmin:dt:Tmax;
N0 = zeros(1,4);
N0(3) = 0.2;
%I =2e-3;TA I =12e-3;
I =0.5e-3;
Vnc= 0.24e-18;
n=3.5;
Nss=(I./(E.*Vp))*tc;
NO = length(omega);
%% main code
%%%%%%%%%%%%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%%%%%%%%%%%
Time = Tmin:dt:Tmax;
z0 = repmat([0.2;0.2;0.2;0.2], 1, NO);
[t,z] = ode45(@rate_eq_program_1,Time,z0);
param_rate_eq
%all the plotting is wrong.
figure(1)
plot(t,z(:,1)); % divided to normalize
xlabel('time [ns]','FontSize',14); % size of x label
ylabel('Arbitrary units','FontSize',14); % size of y label
set(gca,'FontSize',14); % size of tick marks on both axis
legend('N', 'Location','SE') % legend inside the plot
figure(2)
plot(t,z(:,2)); % divided to normalize
xlabel('time [ns]','FontSize',14); % size of x label
ylabel('Arbitrary units','FontSize',14); % size of y label
set(gca,'FontSize',14); % size of tick marks on both axis
legend('NC', 'Location','SE') % legend inside the plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function %%%%%%%%%%%%%%%%%%%%%%%%
function ys = rate_eq_program_1(t,z)
%{
param_rate_eq_fano
%}
z = reshape(z, 4, []);
ys = zeros(size(z));
%rR=(-p*gamma_c)/(1i*(delta_c)+gamma_T-1/2*(1-1i*henry)*conf_NC*V_g*g_n*(z(2)-N_0));
deltawc=KD.*(z(2,:))-1i*KA.*(z(2,:));
deltac=omega0+deltawc-omega;
rR=(-P.*gamac)./(1i*(deltac)+gamat-(1/2*(1-1i)*confinec*vg*gN).*(z(2,:)-Ntr));
rRS=-P.*gamac./(1i.*(deltac)+gamat);
%rRS=rB+(-1i*tB-rB).*gamac./(1i.*(deltac)+gamat);
g=alphai+(1/(2*L)).*log(1./(abs(rRS)).^2);
A=(2*epsilon0*n*ng)./(hbar.*omega0);
B=(1+(abs(rRS))).*(1-(abs(rRS)));
D=((g)-alphai);
H=c./(omega0.*n);
F=imag(rRS)./abs(rRS);
sigmas=A.*((B./D)+(H.*F));
ys(1,:) = I./(E.*Vp)-(z(1,:)./tc)-GN.*(z(1,:)-Ntr).*sigmas.*((abs(z(3,:)).^2)./Vc);
ys(2,:) = (-z(2,:)./tc)-GLA.*(z(2,:)-Ntr).*(abs(z(4,:)).^2)./Vnc;
ys(3,:) = 1/2*(1-1i)*GN.*(z(1,:)-Nss).*z(3,:)+(1/tin.*((z(4,:)./rR)-z(3,:)));
ys(4,:) =((-1i.*deltawc-gamat).*z(4,:))-(P.*gamac.*z(3,:))+(1/2*(1-1i).*confinec*vg*gN.*(z(2,:)-Ntr).*z(4,:));
ys = ys(:);
end
end
  11 Comments
mohammad heydari
mohammad heydari on 4 Feb 2020
thank you very much for your help.
Of course, the disagreement is normal in these cases. And I'll try to include more explanations along with the formulas over the coming hours to give a better description of the situation.
In any case, I sincerely thank you

Sign in to comment.

Categories

Find more on Physics in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by