How to plot a numerical integration?

2 ビュー (過去 30 日間)
Johan Sebastian Diaz Tovar
Johan Sebastian Diaz Tovar 2018 年 8 月 7 日
コメント済み: Aquatris 2018 年 8 月 8 日
Hello guys, I made a numerical integration of this function using adaptive Gauss-Kronrod quadrature 'quadgk()':
as you can see the integrate is k, but for the purpose that I want, I need to plot that result vs z, which is the depth. So, for that I made this code:
%%Total fluence parameters and constants
nu = 0; %zeroth-order Bessel subscript for the function of the first kind.
g = 0.9; %anisotropy coefficient
f = g.^2; g1 = g./(g+1); %first two moments of the delta P1 phase function
n = 1.4; A = -0.13755*n.^3 + 4.3390*n.^2 -4.90366*n + 1.6896; %refraction index and cubic polynomial estimation
%for FRC
Ups = 0.060; Ua = 0.0576; Us = Ups./(1-g); %Reduced, absorption and scattering coefficients
Utr = Ups + Ua; h = 2./(3*Utr);
Ueff = sqrt(3*Ua*Utr); U1s = Us*(1 - f);
U1t = Ua + U1s; z1 = Utr./Ups; %z is the direction of the collimated light within the medium.
r0 = 3./Utr; %Gaussian beam radius.
z = linspace(0,z1,100000);
%%numerical integration
fun = @(x) (((3.*U1s.*(U1t + g.*Ua).*r0.^2.*exp(-(r0.*x).^2)./8)./4.*(x.^2 + Ueff.^2 - U1t.^2)).*exp(-Utr.*z) +
...
((-3.*g1.*U1s.*r0.^2.*exp(-(r0.*x).^2)./8 - 4.*(1./(A.*h) + U1t.^2).* ...
((3.*U1s.*(U1t + g.*Ua).*r0.^2.*exp(-(r0.*x).^2)./8)./4.*(x.^2 + Ueff.^2 - U1t.^2)))./4.*(1./(A.*h) + sqrt(x.^2 +
Ueff.^2))).* ...
exp(-sqrt(x.^2+Ueff.^2).*z)).*besselj(nu,x).*x;
PHId = quadgk(fun,0,Inf);
%%Plot of the total fluence rate, collimated and diffuse fluence rate
PHIc = exp(-U1t.*z);
PHI = PHIc + PHId;
figure;
semilogy(z,PHI)
xlabel('Reduced depth z/l*');
ylabel('Normalized Fluence Rate \phi');
Basically, I made a vector from 0 to z1 (which is the value that I need to reach) and once the numerical integration is made I tried to plot it vs z. I do not know if I am making something wrong, I am new in MATLAB and I am really thankful for your answers.

採用された回答

Aquatris
Aquatris 2018 年 8 月 8 日
The problem is, the quadgk function expect a scalar value from your function. However you assign z to be a vector, hence your function is not scalar. I modified your code and reduce the step size for z to reduce the computational time to be able to test it. Here is my modifications for your code, feel free to ask anything if you do not understand any part.
%%Total fluence parameters and constants
nu = 0; %zeroth-order Bessel subscript for the function of the first kind.
g = 0.9; %anisotropy coefficient
f = g.^2; g1 = g./(g+1); %first two moments of the delta P1 phase function
n = 1.4; A = -0.13755*n.^3 + 4.3390*n.^2 -4.90366*n + 1.6896; %refraction index and cubic polynomial estimation
%for FRC
Ups = 0.060; Ua = 0.0576; Us = Ups./(1-g); %Reduced, absorption and scattering coefficients
Utr = Ups + Ua; h = 2./(3*Utr);
Ueff = sqrt(3*Ua*Utr); U1s = Us*(1 - f);
U1t = Ua + U1s; z1 = Utr./Ups; %z is the direction of the collimated light within the medium.
r0 = 3./Utr; %Gaussian beam radius.
zSet = linspace(0,z1,100);
%%numerical integration
for i = 1:length(zSet)
z = zSet(i);
fun = @(x) (((3.*U1s.*(U1t + g.*Ua).*r0.^2.*exp(-(r0.*x).^2)./8)./4.*(x.^2 + Ueff.^2 - U1t.^2)).*exp(-Utr.*z) + ...
((-3.*g1.*U1s.*r0.^2.*exp(-(r0.*x).^2)./8 - 4.*(1./(A.*h) + U1t.^2).* ...
((3.*U1s.*(U1t + g.*Ua).*r0.^2.*exp(-(r0.*x).^2)./8)./4.*(x.^2 + Ueff.^2 - U1t.^2)))./4.*(1./(A.*h) + sqrt(x.^2 + ...
Ueff.^2))).* ...
exp(-sqrt(x.^2+Ueff.^2).*z)).*besselj(nu,x).*x;
PHId(i) = quadgk(fun,0,Inf);
%%Plot of the total fluence rate, collimated and diffuse fluence rate
PHIc(i) = exp(-U1t.*z);
PHI(i) = PHIc(i) + PHId(i);
end
figure;
semilogy(zSet,PHI)
xlabel('Reduced depth z/l*');
ylabel('Normalized Fluence Rate \phi');
  4 件のコメント
Johan Sebastian Diaz Tovar
Johan Sebastian Diaz Tovar 2018 年 8 月 8 日
Thank you so much!!!! I am really thankful with your answers!
Aquatris
Aquatris 2018 年 8 月 8 日
Happy to help. Hope these are what you were looking for.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeGraphics Performance についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by