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.
0 个评论
采纳的回答
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
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
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
更多回答(1 个)
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
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!