problem with numerical integration

Question

zina shadidi 2021-2-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/751954-problem-with-numerical-integration

编辑： zina shadidi 2021-2-23

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)

2 个评论
显示无隐藏无

Star Strider 2021-2-21

在 MATLAB Online 中打开

Unrecognized function or variable 'fun_disp1'.
  w1 = fun_disp1(k,a) ;

Paul Hoffrichter 2021-2-21

编辑：Paul Hoffrichter 2021-2-21

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 2021-2-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/751954-problem-with-numerical-integration#answer_629509

在 MATLAB Online 中打开

 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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Walter Roberson 2021-2-22

编辑：Walter Roberson 2021-2-22

There are three essential cases for integral():

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)
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)
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 2021-2-23

编辑：zina shadidi 2021-2-23

Thank you very much

very clear explination .

请先登录，再进行评论。

problem with numerical integration

2 个评论
显示无隐藏无

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

problem with numerical integration

2 个评论 显示 无隐藏 无

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

2 个评论
显示无隐藏无

6 个评论
显示 4更早的评论隐藏 4更早的评论