How to integrate a function inside the for loop ?

I am trying to do integration inside a for loop. I am having problem with the matrix dimension or syntax maybe. Would anyone please help me fix it. Also It would like to plot the result of the integration.

clc
close all
clear all

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
Q=1.6*10^-16;                 % Electronic Charge
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy


A=((2*Fs)/(h^3*c^2));
B=(K*Ts);
C=((2*Fa)/(h^3*c^2));
D=(K*Ta);
vdata=0:0.01:1.6;
n=length(vdata);
q=zeros(1,n);

for i=1:n
    V=vdata(i);
    fun=@(E)((q.*E.^2).*((A/(exp(E./B)-1))-(C/(exp(E-(Q.*V))./B)-1)-(C/(exp(E./D)-1))));
    q(i)=integral(fun,Eg,inf);
end 

plot(V,q)
    

Answer 1

[Walter Roberson]

Your E will not be a scalar, so having something "/" a calculation involving E would be invoking the "matrix right divide" operator, not the element-by-element division operator.

However even if you fix that, you have the issue that your function involves multiplying by all of q where q is a vector, so the integral would have to result in a vector (one for each element of q). Because integral() passes in a non-scalar for the parameter, that is incompatible -- unless you use the 'vectorized', true option which would pass in a scalar for E and so allow a vector result the same size as q to be returned. But... your context expects a scalar result.

Each time you include all of q in your fun you are using a q that has at least some of the entries still as their zeros() . Each of those leads to a zero result because the 0 will be multiplying the entire result of the term. But.. that includes the first time, where q is all zero, so the first output would return all 0, and no matter which of those all-zero values you selected to store into q(i), you would be storing a zero, so q would not change. And that implies that all of your results are going to come out zero.

I think you need to reconsider using q inside fun . Perhaps it was intended to be Q ?

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
Q=1.6*10^-16;                 % Electronic Charge
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy


A=((2*Fs)/(h^3*c^2));
B=(K*Ts);
C=((2*Fa)/(h^3*c^2));
D=(K*Ta);
vdata=0:0.01:1.6;
n=length(vdata);
q=zeros(1,n);

for i=1:n
    V=vdata(i);
    fun=@(E)((q.*E.^2).*((A./(exp(E./B)-1))-(C./(exp(E-(Q.*V))./B)-1)-(C./(exp(E./D)-1))));
    whos
    q(i)=integral(fun,Eg,inf);
end 

  Name       Size             Bytes  Class              Attributes

  A          1x1                  8  double                       
  B          1x1                  8  double                       
  C          1x1                  8  double                       
  D          1x1                  8  double                       
  Eg         1x1                  8  double                       
  Fa         1x1                  8  double                       
  Fs         1x1                  8  double                       
  K          1x1                  8  double                       
  Q          1x1                  8  double                       
  Ta         1x1                  8  double                       
  Ts         1x1                  8  double                       
  V          1x1                  8  double                       
  c          1x1                  8  double                       
  fun        1x1                 32  function_handle              
  h          1x1                  8  double                       
  i          1x1                  8  double                       
  n          1x1                  8  double                       
  q          1x161             1288  double                       
  vdata      1x161             1288  double                       

Arrays have incompatible sizes for this operation.

Error in solution (line 22)
fun=@(E)((q.*E.^2).*((A./(exp(E./B)-1))-(C./(exp(E-(Q.*V))./B)-1)-(C./(exp(E./D)-1))));

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

plot(V,q)
    

[Md. Golam Zakaria]

you are right it was supposed to be Q. My typing mistake. But does that solve my problem? I am getting weird graph.

[Walter Roberson]

Your function goes to infinity as E goes to infinity -- at least for positive V it does. So your integrals are always infinite.

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
Q=1.6*10^-16;                 % Electronic Charge
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy


A=((2*Fs)/(h^3*c^2));
B=(K*Ts);
C=((2*Fa)/(h^3*c^2));
D=(K*Ta);
vdata=0:0.01:1.6;
n=length(vdata);
q=zeros(1,n);

syms E V
F = ((Q.*E.^2).*((A./(exp(E./B)-1))-(C./(exp(E-(Q.*V))./B)-1)-(C./(exp(E./D)-1))))

F = 

string(F)

ans = "-(6490371073168535*E^2*((7646125632079076637682471796736*exp((6490371073168535*V)/40564819207303340847894502572032 - E))/15625 - 21313242180011364974592/(exp((15625*E)/7749) - 1) + 986724175000526085647499264/(exp((100000*E)/2583) - 1) - 1))/40564819207303340847894502572032"

L0 = limit(limit(F, V, 0),E,inf)

L0 = 

∞

L1 = limit(limit(F, V, 1),E,inf)

L1 = 

∞

limit(F, E, inf)

ans = 

∞

[Md. Golam Zakaria]

The following code gives a reasonable output. Would you take a look and confirm whether I can use the plot command as I used.

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
Q=1.6*10^-16;                 % Electronic Charge
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy


A=((2*Fs)/(h^3*c^2));
B=(K*Ts);
C=((2*Fa)/(h^3*c^2));
D=(K*Ta);
vdata=0:0.001:1.6;
n=length(vdata);
q=zeros(1,n);

for i=1:n
    V=vdata(i);
    fun=@(E)  Q.*((E.^2).*((A./(exp(E./B)-1))-(C./(exp(E-(V))./B)-1)-(C./(exp(E./D)-1))));
    q(i)=integral(fun,Eg,inf);
end 

plot(vdata,q)
    

[Walter Roberson]

No. Your integral is properly infinite, but integrate() is giving up on integrating it after a while.

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
Q=1.6*10^-16;                 % Electronic Charge
h= 4.136*10^-15;              % Plancks Constant
c= 3*10^8;                    % Speed of light
K = 8.61*10^-5;               % Boltzmanns Constant
Ts=5760;                      % Temparature of the sun
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy


A=((2*Fs)/(h^3*c^2));
B=(K*Ts);
C=((2*Fa)/(h^3*c^2));
D=(K*Ta);
vdata=0:0.001:1.6;
n=length(vdata);
q=zeros(1,n);
syms V E

fun = Q.*((E.^2).*((A./(exp(E./B)-1))-(C./(exp(E-(V))./B)-1)-(C./(exp(E./D)-1))))

fun = 

string(fun)

ans = "-(6490371073168535*E^2*((7646125632079076637682471796736*exp(V - E))/15625 - 21313242180011364974592/(exp((15625*E)/7749) - 1) + 986724175000526085647499264/(exp((100000*E)/2583) - 1) - 1))/40564819207303340847894502572032"

limit(fun, E, inf)

ans = 

∞