How to plot a numerical integration?
3 views (last 30 days)
Show older comments
Johan Sebastian Diaz Tovar
on 7 Aug 2018
Commented: Aquatris
on 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.
0 Comments
Accepted Answer
Aquatris
on 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 Comments
More Answers (0)
See Also
Categories
Find more on Graphics Performance in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!