Issue in sloving numerical integration
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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
figure;
imp1 = plot(rho,phi,'Color','blue','LineStyle','-');
5 Commenti
Torsten
il 25 Nov 2022
Modificato: Torsten
il 25 Nov 2022
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)
It's impossible for us to look through your code and decide which function handle you want to evaluate with which numerical input.
Risposte (0)
Vedere anche
Categorie
Scopri di più su Bessel functions in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!