How to perform integration inside for loop? [matrix dimension must agree issue]

1 次查看(过去 30 天)
Hello Everyone I am a beginner at matlab. I am trying to perform integration inside a for loop. But I am getting 'matrix dimension must agree' error. I have a feeling that this is a minor problem and any expert can solve this problem within minuites. Would anyone be kind enough to spare few minuites to solve this proble. The code is below-
clc
clear all
close all
h=6.582*10^-16;
k=8.617*10^-5;
T=300;
beta=7.021*10^-4;
gamma=1108;
C1=5.5;
C2=4;
A1=3.231*10^2;
A2=7.237*10^3;
Ad=1.052*10^6;
Ep1=1.82*10^-2;
Ep2=5.773*10^-2;
Egd=3.2;
Eg0_1=1.1557;
Eg0_2=2.5;
Eg1=Eg0_1-((beta*(T^2))/(T+gamma));
Eg2=Eg0_2-((beta*(T^2))/(T+gamma));
Fs= 2.16*10^-5*pi; % Geometrical Factor for sun
H= 4.136*10^-15; % Plancks Constant
c= 3*10^8; % Speed of light
K = 8.6173*10^-5; % Boltzmanns Constant
Ts=5760; % Temparature of the sun
q=1.6*10^-19;
A=((2*Fs)/((H^3)*(c^2)));
x=1:0.004:5;
num=numel(x);
output=nan(1,num);
for v=1:num
lambda=0.2:0.001:1.2;
Irradiance=(A.*(((H*c)./lambda).^3./(exp((((H*c)./lambda)./(K.*Ts)))-1))).*q;
alpha=C1*A1*(((((h.*((2*pi)./lambda))-Eg1+Ep1).^2)./(exp(Ep1/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg1-Ep1).^2)./(1-exp(-Ep1/(k*T)))))+C2*A2*(((((h.*((2*pi)./lambda))-Eg2+Ep2).^2)./(exp(Ep2/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg2-Ep2).^2)./(1-exp(-Ep2/(k*T)))))+Ad.*((2*pi.*(c./lambda))-Egd).^(1/2);
depth=v;
attenuation=@(depth) (depth.*alpha);
alpha_x=integral(attenuation,0,depth');
output(v)=alpha.*Irradiance.*exp(-alpha_x);
end
plot(x,output)
For further simplification, I should note that I triend to use nested for loop for lambda, but that complecated the problem for me. But this code should also work.

采纳的回答

VBBV
VBBV 2022-2-26
编辑:VBBV 2022-2-26
clc
clear all
close all
h=6.582*10^-16;
k=8.617*10^-5;
T=300;
beta=7.021*10^-4;
gamma=1108;
C1=5.5;
C2=4;
A1=3.231*10^2;
A2=7.237*10^3;
Ad=1.052*10^6;
Ep1=1.82*10^-2;
Ep2=5.773*10^-2;
Egd=3.2;
Eg0_1=1.1557;
Eg0_2=2.5;
Eg1=Eg0_1-((beta*(T^2))/(T+gamma));
Eg2=Eg0_2-((beta*(T^2))/(T+gamma));
Fs= 2.16*10^-5*pi; % Geometrical Factor for sun
H= 4.136*10^-15; % Plancks Constant
c= 3*10^8; % Speed of light
K = 8.6173*10^-5; % Boltzmanns Constant
Ts=5760; % Temparature of the sun
q=1.6*10^-19;
A=((2*Fs)/((H^3)*(c^2)));
x=1:0.04:5;
num=numel(x);
output=zeros(1,num);
lambda=0.2:0.01:1.2
lambda = 1×101
0.2000 0.2100 0.2200 0.2300 0.2400 0.2500 0.2600 0.2700 0.2800 0.2900 0.3000 0.3100 0.3200 0.3300 0.3400 0.3500 0.3600 0.3700 0.3800 0.3900 0.4000 0.4100 0.4200 0.4300 0.4400 0.4500 0.4600 0.4700 0.4800 0.4900
1:1.2;
for v=1:num
Irradiance=(A.*(((H*c)./lambda).^3./(exp((((H*c)./lambda)./(K.*Ts)))-1))).*q;
alpha=C1*A1*(((((h.*((2*pi)./lambda))-Eg1+Ep1).^2)./(exp(Ep1/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg1-Ep1).^2)./(1-exp(-Ep1/(k*T)))))+C2*A2*(((((h.*((2*pi)./lambda))-Eg2+Ep2).^2)./(exp(Ep2/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg2-Ep2).^2)./(1-exp(-Ep2/(k*T)))))+Ad.*((2*pi.*(c./lambda))-Egd).^(1/2);
attenuation=@(depth) (depth.*alpha);
depth=v;
alpha_x(v,:)=integral(@(depth) attenuation(depth),0,depth,'ArrayValued',true);
output(v,:)=alpha_x(v,:).*Irradiance.*exp(-depth);
end
plot(x,output)
axis([1 2 0 2600])
  2 个评论
VBBV
VBBV 2022-2-26
编辑:VBBV 2022-2-26
When computing output , why do you use alpha instead of alpha_x ?
What is the use of evaluating alpha_x that is not used anywhere in the code ?
Put
depth = v
after
attenuation=@(depth) (depth.*alpha);
This allows the assigned value for depth to be used correctly, otherwise it uses only symbolic variable
Md. Golam Zakaria
Md. Golam Zakaria 2022-2-26
@VBBV actually the first code I used in the question thread had a mistake.The right code is
clc
clear all
close all
h=6.582*10^-16;
k=8.617*10^-5;
T=300;
beta=7.021*10^-4;
gamma=1108;
C1=5.5;
C2=4;
A1=3.231*10^2;
A2=7.237*10^3;
Ad=1.052*10^6;
Ep1=1.82*10^-2;
Ep2=5.773*10^-2;
Egd=3.2;
Eg0_1=1.1557;
Eg0_2=2.5;
Eg1=Eg0_1-((beta*(T^2))/(T+gamma));
Eg2=Eg0_2-((beta*(T^2))/(T+gamma));
Fs= 2.16*10^-5*pi; % Geometrical Factor for sun
H= 4.136*10^-15; % Plancks Constant
c= 3*10^8; % Speed of light
K = 8.6173*10^-5; % Boltzmanns Constant
Ts=5760; % Temparature of the sun
q=1.6*10^-19;
A=((2*Fs)/((H^3)*(c^2)));
x=1:0.004:5;
num=numel(x);
output=nan(1,num);
for v=1:num
lambda=0.2:0.001:1.2;
Irradiance=(A.*(((H*c)./lambda).^3./(exp((((H*c)./lambda)./(K.*Ts)))-1))).*q;
alpha=C1*A1*(((((h.*((2*pi)./lambda))-Eg1+Ep1).^2)./(exp(Ep1/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg1-Ep1).^2)./(1-exp(-Ep1/(k*T)))))+C2*A2*(((((h.*((2*pi)./lambda))-Eg2+Ep2).^2)./(exp(Ep2/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg2-Ep2).^2)./(1-exp(-Ep2/(k*T)))))+Ad.*((2*pi.*(c./lambda))-Egd).^(1/2);
depth=v;
attenuation=@(depth) (depth.*alpha);
alpha_x=integral(attenuation,0,depth');
output(v)=alpha.*Irradiance.*exp(-alpha_x);
end
plot(x,output)
@VBBV please reflect your advice on this one. As you mentioned I was not using alpha_x . Its was supposed to be inside the exp. Would you please redo the whole thing again?please

请先登录,再进行评论。

更多回答(1 个)

Torsten
Torsten 2022-2-25
Use
alpha_x{v} = integral(attenuation,0,depth,'ArrayValued',true);
output{v} = alpha.*Irradiance.*exp(-depth);
instead of
alpha_x=integral(attenuation,0,depth');
output(v)=alpha.*Irradiance.*exp(-depth);
But the result of your integration will simply be
alpha_x{v} = v*alpha
Is this really what you want ?
  3 个评论
Torsten
Torsten 2022-2-25
编辑:Torsten 2022-2-25
Yes, you need curly brackets for alpha_x as well as for output.
The reason is that for each v, alpha_x is a vector of the same length as lambda. In order to save these "num" vectors, they have either to be saved in a cell array (as done above) or in a matrix
alpha_x(:,v) = integral(attenuation,0,depth,'ArrayValued',true);
output(:,v) = alpha.*Irradiance.*exp(-alpha_x);
But the main question is whether the integral really gives the answer you expect, since - as said above -
alpha_x{v} = v*alpha
So no integration is needed to get this result.
Md. Golam Zakaria
Md. Golam Zakaria 2022-2-25
编辑:Md. Golam Zakaria 2022-2-25
@Torsten Actually some books use integration some dont.I skipped the R(E) part.
clc
clear all
close all
h=6.582*10^-16;
k=8.617*10^-5;
T=300;
beta=7.021*10^-4;
gamma=1108;
C1=5.5;
C2=4;
A1=3.231*10^2;
A2=7.237*10^3;
Ad=1.052*10^6;
Ep1=1.82*10^-2;
Ep2=5.773*10^-2;
Egd=3.2;
Eg0_1=1.1557;
Eg0_2=2.5;
Eg1=Eg0_1-((beta*(T^2))/(T+gamma));
Eg2=Eg0_2-((beta*(T^2))/(T+gamma));
Fs= 2.16*10^-5*pi; % Geometrical Factor for sun
H= 4.136*10^-15; % Plancks Constant
c= 3*10^8; % Speed of light
K = 8.6173*10^-5; % Boltzmanns Constant
Ts=5760; % Temparature of the sun
q=1.6*10^-19;
A=((2*Fs)/((H^3)*(c^2)));
x=5;
lambda=0.2:0.001:1.2;
alpha=C1*A1*(((((h.*((2*pi)./lambda))-Eg1+Ep1).^2)./(exp(Ep1/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg1-Ep1).^2)./(1-exp(-Ep1/(k*T)))))+C2*A2*(((((h.*((2*pi)./lambda))-Eg2+Ep2).^2)./(exp(Ep2/(k*T))-1))+((((h.*((2*pi)./lambda))-Eg2-Ep2).^2)./(1-exp(-Ep2/(k*T)))))+Ad.*((2*pi.*(c./lambda))-Egd).^(1/2);
Irradiance=(A.*(((H*c)./lambda).^3./(exp((((H*c)./lambda)./(K.*Ts)))-1))).*q;
G=alpha.*Irradiance;
plot(lambda,G)
For example in this code if I skip the exponential part I get a qualititively accurate graph. But when I add the exponential part I get zero. Thats why I thought maybe I need to integrate to get non zero value.
can you please compile the full code with the suggestion you made and see If you get an output graph.

请先登录,再进行评论。

类别

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

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by