Issue in sloving numerical integration
조회 수: 1 (최근 30일)
이전 댓글 표시
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
Torsten
2022년 11월 25일
편집: Torsten
2022년 11월 25일
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.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Bessel functions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!