How to perform integration inside for loop? [matrix dimension must agree issue]
Show older comments
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.
Accepted Answer
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
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 Comments
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
on 26 Feb 2022
More Answers (1)
Torsten
on 25 Feb 2022
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 Comments
Md. Golam Zakaria
on 25 Feb 2022
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
on 25 Feb 2022
Edited: Md. Golam Zakaria
on 25 Feb 2022
@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.
Categories
Find more on Mathematics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
