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

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.

Answer 1

[VBBV]

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

[VBBV]

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]

@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

Answer 2

[Torsten]

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 ?

[Md. Golam Zakaria]

@Torsten actually I had aslight mistake in my code. The code was supposed to be

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)

now as you have answerd , do I have to use curly brackets?

[Torsten]

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]

@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.