How to plot a numerical integration?

2 Ansichten (letzte 30 Tage)
Johan Sebastian Diaz Tovar
Kommentiert: Aquatris am 8 Aug. 2018
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.

Akzeptierte Antwort

Aquatris
Aquatris am 8 Aug. 2018
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 Kommentare
Johan Sebastian Diaz Tovar
Thank you so much!!!! I am really thankful with your answers!
Aquatris
Aquatris am 8 Aug. 2018
Happy to help. Hope these are what you were looking for.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Graphics Performance finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by