problem with numerical integration

I have problem with integration in the folowing code at the last line

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phonon dispersion relation and density of state for a simple monoclinic
 % lattice using the lenier spring model
 % input parameters 
 close all; clear all;
 k_B = 1.38e-23; % Boltzmann's constant
h_bar = 6.626e-34/2/pi; % Reduced Planck's constant
T=300;
 d = 3 ; % dimensions in NiTi 
m2=47.867;  % T1 atomic mass in U
m1=58.693; % Ni atomic mass in U
 a =2.89 ; %nearest neighbour in real space monoclinic NiTi

 % quotient of the spring constamts for first &second nearest neibours
 % c-quot = c1/c2
c_quot = [inf , 6 , 3 , 2 , 1 ] ; %( guess )

 % number of point sink ?space
 n = 1000 ; % ploting dispersion relation
 n_dos = 10000000 ; % calculating  DoS

 % number of histogram bins for DoS
w_bin =150 ;

 % brilion zone symmetry point
 Z = [0;0;0] ;
 D=[1/2;-1/2;0];
Y = [0;1/2;0] ;
G=[0;0;0];
 X = [1/2 ; 0 ; 0 ] ;
 A = [ 0 ; 0 ; 1/2] ;

    
 %dispersion relation plot order (symmetry points)
 po= [ Z,D,Y,G ,X,A] ;
 po_label = {'Z','D','Y','G','X','A'};
 
 %phonon dispersion relation
 %creating n linearly spaced k vector between each pear of symetry
 n_po = size(po,2)-1; % M = size(X,DIM) returns the length of the dimension specified
    %by the scalar DIM.
 k = nan(d,(n-1)*n_po+1) ; % kx ; ky ; kz (d by (n-1)*n-p0+1)) matrix of undefined  values
 kk = zeros(1 ,size(k ,2 ));  % length of path in k- space for plot axis
 kk_s= zeros ( 1, n_po+1) ; % value of kk at the symetry points
 %zeros is a matrix all its elements equal zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 k(:,1)= po(:,1) ;
 kk(1) = kk_s(1) ;
 for ind1=1:n_po % loop from 1 to the matrix size

 ind1_start =( ind1-1)*( n-1)+2; % dispertion relation equations
 ind1_end =ind1*( n-1)+1;        % dispertion relation equations

 kk_s(ind1 +1)= kk_s( ind1)+norm( po (:,ind1+1)-po( :,ind1)) ;
 temp_kk = linspace(kk_s(ind1),kk_s(ind1 +1),n) ;
 %linspace(X1, X2,n) generates a row vector of n linearly
  %equally spaced points between X1 and X2.
 
 kk(1,ind1_start :ind1_end)=temp_kk(2:end);

 for ind2=1:d
 temp_k = linspace(po(ind2,ind1 ),po(ind2,ind1+1),n) ;
 k(ind2,ind1_start:ind1_end)= temp_k(2:end);
 end
end

 % calculating the frequencies  omega /(sqrt(C/m)) for each k  from the dispersion relation:
 % nearest neigbhours
 
  w1 = fun_disp1(k,a) ;
 
 if k == 0
         w1= 0;
elseif k == pi
         w1 =2* abs(w1/sqr(m1));
end
 
 
 % nearest & next nearest neighbours
 numel_c_quot = numel(c_quot);  %N = numel(A) returns the number of elements, N, in array A
 w2 = nan( size(k ,1) ,size(k,2),numel_c_quot ) ;
 for ind1 =1:numel_c_quot
 w2(:,:,ind1)=fun_disp2(k,a,c_quot(ind1));
 if w2 == 0;
         w2= 2*w2*abs((sqr(1/m1+1/m2)));
       elseif w2 == pi
        w2 = 2*w2*abs((sqr(1/m2)));
 end
 end



 %ploting the dispersion relation 
 % nearest neighbours
 ymax1 = round(11*max(w1(:)))/10;

figure
 plot(kk,w1,'b-' )
for ind1= 1:(n_po+1)
 line([kk_s(ind1) kk_s(ind1)],[0 ymax1 ] , 'Color' ,' r ' )
 end
 ylim([0 ymax1])
 xlim ([ kk(1)  kk(end-1)])
 %'Phonon dispersion relation, nearest neighbours=PDR 
 % nearest neighbours=1st
  title( 'Phonon dispersion relation of monoclinic NiTi , nearest neighbours' )

 n=13; ylabel([ ' frequency^{' num2str(n) '}','Hz' ]);
  xlabel ( ' (wavevector) ' )
 set (gca , 'xtick' , kk_s )  
 %gca returns the current axes or chart for the current figure (plot information in xlim,ylim,title)
 xticklabels({'Z','D','Y','G','X','A'})
  set (gca   ,'xticklabels', xticklabels )

 %nearest & next nearest neighbours
 for ind1 = 1 : numel_c_quot
 ymax2 = round(11*max(max(w2)))/10;

figure
 plot(kk,w2(:,:,ind1),'b-')
 for ind2 = 1:( n_po+1)
 line([kk_s(ind2) kk_s( ind2)] , [0 ymax2( ind1)] )
 end
 ylim([0 ymax2(ind1)])
 xlim([kk(1) kk(end)])
  title( 'Phonon dispersion relation of monoclinic NiTi , 1st &2nd nearest neighbours' )
 n=13; ylabel([ ' frequency^{' num2str(n) '}','Hz' ]);
  xlabel ( ' (wavevector) ' )
 
 end


%density of state-------


% choosing random K vectors in the first brillion zone
k_rand = 2*pi /a*( rand(d ,n_dos )-0.5) ;
% calculating the frequences from the dispersion relation
% nearest neighbours
w1_rand = fun_disp1 (k_rand , a ) ;
if k == 0
         w1_rand = 0;
elseif k == pi
         w1 = 2*w1/abs(sqr(m1));
end

% nearest &next nearest neighbours
 w2_rand = nan( size( k_rand ,1 ) ,size( k_rand ,2),numel_c_quot ) ;
 if w2_rand == 0;
         w2_rand =2* w2_rand*(abs(sqr(1/m1+1/m2)));
       elseif w2_rand == pi
        w2_rand = 2*w2_rand(abs(sqr(1/m2)));
     
end 

 for ind1 = 1 : numel_c_quot
 w2_rand ( : , : , ind1 ) = fun_disp2 ( k_rand , a , c_quot ( ind1 ) ) ;
 end
% histogram
[ dw1 ,wn1 ] = hist( w1_rand( : ), w_bin ) ;
% N = hist(Y) bins the elements of Y into 10 equally spaced containers
    %and returns the number of elements in each container.  If Y is a
   % matrix, hist works down the columns.[N,X] = hist(...) also returns the position of the bin centers in X.
 
  % hist(...) without output arguments produces a histogram bar plot of the results.
  %The bar edges on the first and last bins may extend to cover the min and max of the data unless a matrix of data is supplied.
 dos1 = dw1/ n_dos* w_bin /(max(wn1)-min(wn1) )* 1/( a ^3) ; % normalisation
 
 
 % % nearest & next nearest neighbours
 dw2 = nan ( numel_c_quot , w_bin ) ;
 wn2 = nan ( numel_c_quot , w_bin ) ;
 dos2 = nan ( numel_c_quot , w_bin ) ;
 for ind1 = 1 : numel_c_quot
w2_rand_ind1 = w2_rand ( : , : , ind1 ) ;
[ dw2( ind1 , : ) ,wn2( ind1 , : ) ] = hist ( w2_rand_ind1( : ) , w_bin ) ;
dos2( ind1,:)=dw2(ind1,:)/n_dos-w_bin/(max(wn2(ind1,:))-min(wn2(ind1,:)))*1/( a ^3) ;
 %normalisation
k_B = 1.38e-23; % Boltzmann's constant
h_bar = 6.626e-34/2/pi; % Reduced Planck's constant
cv=(1/4*k_B*T)*((h_bar*wn2(ind1,:)).^2/(sinh(h_bar*wn2(ind1,:)/2*k_B*T)).^2).*dos2(ind1,:);

% Calculate the low temperature approx. of the Debye model:

 TD=461.75;
 omega_D = k_B * TD / h_bar;
 
 fun=@(w2)(dos2(ind1,:).*log(2*(sinh(h_bar*wn2(ind1,:)./2*k_B*T))));

 E=(n_dos*k_B*T).*integral(fun,max(wn2(ind1,:)),min(wn2(ind1,:)))
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

the error is

%%%%%%%%%%%%%%

Error in integralCalc/iterateScalarValued (line 323)

fx = FUN(t).*w;

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)

Q = integralCalc(fun,a,b,opstruct);

Error in cvtry (line 192)

E=integral(fun, 0, omega_D)

[Star Strider]

Unrecognized function or variable 'fun_disp1'.

  w1 = fun_disp1(k,a) ;

[Paul Hoffrichter]

You should put your code into a code block so that it is easy to copy, so that it looks like MATLAB code, and so that the line numbers might match the line numbers in your error message. In the CODE section, highlight your code, and hit the "Insert a line of code" button.

Answer 1

[Walter Roberson]

 fun=@(w2)(dos2(ind1,:).*log(2*(sinh(h_bar*wn2(ind1,:)./2*k_B*T))));
 E=integral(fun, 0, omega_D)

integral() calls the given function handle passing in a variable number of values. With the default parameters to integral(), the output must be something that is the same size as the input; the computations must be independent.

However, your fun takes that variable amount of data... ignores it, and always returns the same size of value, whatever that works out to be.

[zina shadidi]

Thank you Walter Roberson . Still I have ,any help please.

the last steps changed to the folowing

k_B = 1.38e-23; % Boltzmann's constant
h_bar = 6.626e-34/2/pi; % Reduced Planck's constant
cv=(1/4*k_B*T)*((h_bar*wn2(ind1,:)).^2/(sinh(h_bar*wn2(ind1,:)/2*k_B*T)).^2).*dos2(ind1,:);
dT=3;

 TD=461.75;
 Tv = dT:dT:600; 
  fun=@(x)(dos2(ind1,:).*log(2*(sinh(h_bar*wn2(ind1,:)./2*k_B*T))));

omega_D = k_B * TD / h_bar;
E = zeros(size(Tv));

for j = 1:length(Tv)
T = Tv(j);
x_D = h_bar*omega_D/k_B*T;
 
 E(j)=(n_dos*k_B*T).*integral(fun,0,x_D)

end
 end

[zina shadidi]

the error is

Matrix dimensions must agree.

Error in integralCalc/iterateScalarValued (line 323)

fx = FUN(t).*w;

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)

Q = integralCalc(fun,a,b,opstruct);

Error in cvtry (line 195)

E(j)=(n_dos*k_B*T).*integral(fun,0,x_D)

[Walter Roberson]

fun=@(x)(dos2(ind1,:).*log(2*(sinh(h_bar*wn2(ind1,:)./2*k_B*T))));

Where does the calculation use the value x that was passed in? Everything in the body of the anonymous function is already calculated, effectively a constant. You are doing the equivalent of

FUN = (dos2(ind1,:).*log(2*(sinh(h_bar*wn2(ind1,:)./2*k_B*T))));
fun = @(x) FUN;

if that is deliberate, if you are deliberately calculating a scalar constant and want to integrate the constant then you would need.

fun = @(x) FUN.*ones(size(x))

but in such a case why not just say that the integral is (upper minus lower) times FUN??

[zina shadidi]

i want to integrate the folowing function over

 (dos2(ind1,:)

or over T.

the function

FUN = (dos2(ind1,:).*log(2*(sinh(h_bar*wn2(ind1,:)./2*k_B*T))))

the limits is from zero to the final value of

(dos2(ind1,:)

in the case of the variable T I have the integration limits according to the folowing

 T=dT:dT:600

, where

dT=4

and the (dos2(ind1,:) size is (5x150)

[Walter Roberson]

There are three essential cases for integral():

1. You are integrating a constant, such as 0. In such a case, you need to arrange that the output is the same size as the input. For example integral(@(t) 5*ones(size(t)), lower_bound, upper_bound)
2. You are integrating an expression in a single variable and if the expression is passed a scalar then it would return a scalar. . In such a case, the function must be vectorized and return a result the same size as its input. For example integral(@(t) sin(exp(t)), lower_bound, upper_bound) . The expression should involve the variable named in the @() portion of the function handle or the value of the first parameter passed in if you do not use an anonymous function such as integral(@sin, lower_bound, upper_bound)
3. You are integrating an expression in a single variable and for any given value a non-scalar is returned . In such a case, the function does not need to be vectorized, but the 'ArrayValued', true option must be used. The expression should involve the variable named in hte @() portion of the function handle or the value of the first parameter passed in if you do not use an anonymous function. For example integral(@(t) [cos(t), cos(2*t)], lower_bound, upper_bound, 'arrayvalued', true) .

What you have coded does not follow any of those. Your expression does not use the variable named in the @() anonymous function, and your expression does not return a constant the same size as the input, and your expression does not return a scalar, and your expression does not use 'arrayvalued' to indicate that you have more than one output per input.

[zina shadidi]

Thank you very much

very clear explination .