Issue in sloving numerical integration

Question

Captain Rituraj Singh 2022-11-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1862118-issue-in-sloving-numerical-integration

编辑： Torsten 2022-11-25

Running this code I get "ImV = function_handle with value: @(r)-ft1-ft2" instead of number. I want to plot ImV Vs rho but I am not getting the numerical values. I don't know where I am making mistake, please help out. Any help shall be highly appreciated. Thanks in advance

clear all; close all; clc
    
    fmGeV_scale = 5.0574;
    sigma = 0.15; 			
    Nc = 3;
    Nf = 3;
    x = 0.15:0.1:0.5;
    s = size(x,2);
    
    rho = 0.3:0.1:2.0;
    
    for i = 1:s
        T = x(i);
        Lambda_MS = 0.176; %GeV
        %alpha = 0.48;
        
        al1 = 33 - 2*Nf;
        al2 = log(2*pi*T/Lambda_MS);
        al3 = al1*al2;
        alpha = 6*pi/al3;
        
        m_pi2 = 0.0196; %in GeV^2 and 19600 in MeV^2
        
        b1 = 5*m_pi2;% When B = 5mpi2
        b2 = 15*m_pi2;% When B = 15mpi2
        b3 = 25*m_pi2;% When B = 25mpi2
        
        n1 = Nc/3;
        n2 = Nf/6;
        n3 = 4*pi*alpha;
        n4 = n3*(n1 + n2);
        
        d1 = 3*b1*alpha/(2*pi*T.^2);% When B = 5mpi2
        d2 = 3*b2*alpha/(2*pi*T.^2);% When B = 15mpi2
        d3 = 3*b3*alpha/(2*pi*T.^2);% When B = 25mpi2
        
        
        md10(:,i) = T.*sqrt(n4); % Our old form of Debye Mass, B = 0
        md11(:,i) = T.*sqrt(n4 + d1);% When B = 5mpi2
        md12(:,i) = T.*sqrt(n4 + d2);% When B = 15mpi2
        md13(:,i) = T.*sqrt(n4 + d3);% When B = 25mpi2
        
        
        n5 = n3.*n1;
        md20(:,i) = T.*sqrt(n5); % New form of Debye mass, B = 0
        md21(:,i) = T.*sqrt(n5 + d1);% When B = 5mpi2
        md22(:,i) = T.*sqrt(n5 + d2);% When B = 15mpi2
        md23(:,i) = T.*sqrt(n5 + d3);% When B = 25mpi2
        
        
       if T == x(2)
           
           TT = T
           
        r = rho*fmGeV_scale;
             rr = rho;
             
             %first term of the potential
             v20 = alpha*exp(-md20(:,i).*r)./r;% When B = 0
             v21 = alpha*exp(-md21(:,i).*r)./r;% When B = 5mpi2
             v22 = alpha*exp(-md22(:,i).*r)./r;% When B = 15mpi2
             v23 = alpha*exp(-md23(:,i).*r)./r;% When B = 25mpi2
             
    
      
             %second term of the potential
             dc20 = alpha*md20(:,i);% When B = 0
             dc21 = alpha*md21(:,i);% When B = 5mpi2
             dc22 = alpha*md22(:,i);% When B = 15mpi2
             dc23 = alpha*md23(:,i);% When B = 25mpi2
             
    
             
             %third term of the potential
             g2 = 4*pi*alpha;
             mdg2 = 0.33*g2*Nc*(T^2); 
             mu0 = mdg2*sigma/alpha;
             mu = mu0^(1/4);
             var = sqrt(2).*mu.*r;
             bes = besselk(-1/2,var); % Bessel Function
             gm1 = gamma(1/4); %Gamma(1\4)
             nom = gm1*sigma*bes;%nominator
             dnom = sqrt(pi)*mu*(2^(3/4));%denominator
             v30 = nom/dnom;
             
    
             
    
             gm2 = gamma(3/4); %Gamma(3\4)
             v40 = gm1*sigma/(2*gm2*mu);
             
     
     
     syms x5 %calculation for g(mDr)
    
      h1 = @(x5) (x5./(x5.*x5 + 1));
      h2 = @(x5,r) (1 - sin(md21.*r.*x5)./(md21.*r.*x5));
      h   = @(x5,r) h1(x5) .* h2(x5,r);
    
      g_mDr = @(r) integral(@(x5) h(x5), 0, inf);
      
      
     
     
     syms x1 r1 %calculation for g(mDx)
    
      h1 = @(x1) (x1./(x1.*x1 + 1));
      h2 = @(x1,r1) (1 - sin(md21.*r1.*x1)./(md21.*r1.*x1));
      h   = @(x1,r1) h1(x1) .* h2(x1,r1);
    
      g_mD = @(r1) integral(@(x1) h(x1,r1), 0, inf);
      
      mr = sqrt(2)*mu.*r;
      Dj1 = besselk(-1/2,mr);
      im = 1i;
      
      Dj2 = real(besselk(-1/2,im*mr));% real part of imaginary bessel function
      
      Dj3 = besselk(-1/2,0);
      
      
     syms x2
     
       mmu = sqrt(2).*mu
       Dx1 =@(x2) besselk(-1/2,mmu.*x2);
       Dx2 =@(x2) real(besselk(-1/2,im*mmu.*x2));
    
    % first integrand and integral
    
    Im1 = @(x2) Dx2(x2).*x2.*x2.*g_mD(x2);
    ph1 =@(r) Dj1.*integral(Im1(x2),0,r);
    
    Im2 = @(x2) Dx1(x2).*x2.*x2.*g_mD(x2);
    ph2 =@(r) Dj2.*integral(Im1(x2),r,inf);
     
    
    Im3 = @(x2) Dx1(x2).*x2.*x2.*g_mD(x2);
    ph3 =@(r) Dj3.*integral(Im1(x2),0,inf);
    
    % r = 0.1:0.1:2.0;
    
    phi =@(r) ph1(r) + ph2(r) - ph3(r);
    
    
    ft1 =@(r) 2*alpha*T*g_mDr;
    ft2 =@(r) mdg2*sigma*T*phi/mu;
    
    ImV =@(r) -ft1-ft2
             
       end
        
    end   
TT = 0.2500
mmu = 0.8265
ImV = function_handle with value:
    @(r)-ft1-ft2
        
    
    
     figure;
    
     imp1 = plot(rho,phi,'Color','blue','LineStyle','-');
Error using plot
Invalid data argument.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Captain Rituraj Singh 2022-11-25

在 MATLAB Online 中打开

Even after making the suggested changes I am still getting same output:

ImV = function_handle with value: @(r)-ft1(r)-ft2(r)

while I need, proper solution in numbers not in algebric format, please help me in that.

Thanks

Torsten 2022-11-25

编辑：Torsten 2022-11-25

在 MATLAB Online 中打开

A function handle is a function. It is necessary to give numerial input to a function handle to get numerical output.

f = @(x) x.^2;
x = 0:0.1:1;
f(x)
ans = 1×11
         0    0.0100    0.0400    0.0900    0.1600    0.2500    0.3600    0.4900    0.6400    0.8100    1.0000

It's impossible for us to look through your code and decide which function handle you want to evaluate with which numerical input.

请先登录，再进行评论。

请先登录，再回答此问题。