How to plot a numerical integration?

Question

Johan Sebastian Diaz Tovar 2018-8-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/413906-how-to-plot-a-numerical-integration

评论： Aquatris 2018-8-8

采纳的回答： Aquatris

在 MATLAB Online 中打开

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Aquatris 2018-8-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/413906-how-to-plot-a-numerical-integration#answer_331970

在 MATLAB Online 中打开

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');