Unable to perform assignment because the left and right sides have a different number of elements in the ode45
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
mohammad heydari
el 1 de Feb. de 2020
Comentada: mohammad heydari
el 4 de Feb. de 2020
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 comentarios
Walter Roberson
el 1 de Feb. de 2020
You should format your code as a block to make it easier for people to copy it.
Walter Roberson
el 1 de Feb. de 2020
What is param_rate_eq_fano ? Is it a script that defines the n used in
A=(2*epsilon0*n*ng)./(hbar.*omega0);
Respuesta aceptada
Walter Roberson
el 1 de Feb. de 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:
- 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
- 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 comentarios
Más respuestas (1)
Walter Roberson
el 2 de Feb. de 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 comentarios
Walter Roberson
el 4 de Feb. de 2020
Sorry, I do not appear to be able to communicate successfully with you on this matter. Perhaps some other volunteer will have more success. However, in my experience, most of the other volunteers are much stricter about requiring commented code than I am.
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!