Can't integrate function using Matlab

I am trying to do a basic integration of a formula for 4 different values. After doing a bit of reading, apparently I need to create a function handle but I am not getting any success. I have put my attempt below. I might be missing more though.

L=0:100;
L_bar=25;
sigma_d =[0.25, 0.5, 0.75, 1.0]; 
mu = log(L_bar) - 0.5*(sigma_d).^2;
L_tilda = exp(mu);

% eta_r = P.*exp(-r./L);
%  I want to integrate eta from zero to infinity where P has changing
%  variables to give 4 different plots.


r = 0:100;

% I was just testing below to see if a loop would run before attempting to
% plot anything

P = zeros(1,length(L));
fun = @(r) P.*exp(-r./L);
for ii=1:numel(sigma_d)
P= 1./(L.*sqrt(2*pi)*sigma_d(ii)).*exp(-log(L/L_tilda(ii)).^2/(2*sigma_d(ii)^2)) ; 
eta_r = integral(fun,0,inf);
end 

Answer 1

[Torsten]

L_bar=25;
sigma_d =[0.25, 0.5, 0.75, 1.0]; 
mu = log(L_bar) - 0.5*(sigma_d).^2;
L_tilda = exp(mu);
r = 0:100;
for ii=1:numel(sigma_d)
    P= @(L)1./(L.*sqrt(2*pi)*sigma_d(ii)).*exp(-log(L/L_tilda(ii)).^2/(2*sigma_d(ii)^2));
    fun = @(L,r) P(L).*exp(-r./L);
    eta_r(ii,:)=integral(@(L)fun(L,r),0,Inf,'ArrayValued',true);
end
plot(r,eta_r)
grid on

[David Harra]

This gives me the exact plots that I was looking for. Perfect

[David Harra]

Torsten I have one question about the plot - only the appearance . I am trying ti change each line with a different line style. For some reason all the lines are plotting the same and I have defined them to be different. using the variable LineTypes

L=0:100;
L_bar=25;
sigma_d =[0.25, 0.5, 0.75, 1.0]; 
mu = log(L_bar) - 0.5*(sigma_d).^2;
L_tilda = exp(mu);
r = 0:100;
%  LineTypes={'-'	'--' ':' '-.'};
for ii=1:numel(sigma_d)
    P= @(L)1./(L.*sqrt(2*pi)*sigma_d(ii)).*exp(-log(L/L_tilda(ii)).^2/(2*sigma_d(ii)^2));
    fun = @(L,r) P(L).*exp(-r./L);
    eta_r(ii,:)=integral(@(L)fun(L,r),0,Inf,'ArrayValued',true);
end
figure(4)
LineTypes={'-'	'--' ':' '-.'};
plot(r,eta_r,'LineWidth',2 , 'LineStyle',LineTypes{ii})
xlabel('r (\mu m)')
xticks([0 20 40 60 80 100])
ylabel('\eta(r)');
ylim([0 1])
title('Spatial Correlation Function LBAR=25')
set(gca, 'FontSize',10);
legend('\sigma_d =0.25', '\sigma_d= 0.5', '\sigma_d =0.75', '\sigma_d = 1')

[Dyuman Joshi]

You will have to apply that individually -

L=0:100;
L_bar=25;
sigma_d =[0.25, 0.5, 0.75, 1.0]; 
mu = log(L_bar) - 0.5*(sigma_d).^2;
L_tilda = exp(mu);
r = 0:100;
LineTypes={'-'	'--' ':' '-.'};
for ii=1:numel(sigma_d)
    P= @(L)1./(L.*sqrt(2*pi)*sigma_d(ii)).*exp(-log(L/L_tilda(ii)).^2/(2*sigma_d(ii)^2));
    fun = @(L) P(L).*exp(-r./L);
    eta_r(ii,:)=integral(@(L)fun(L),0,Inf,'ArrayValued',true);
    plot(r,eta_r(ii,:),'LineStyle', LineTypes{ii}, 'LineWidth', 1.5)
    hold on
end
xlabel('r (\mu m)')
xticks([0 20 40 60 80 100])
ylabel('\eta(r)');
ylim([0 1])
title('Spatial Correlation Function LBAR=25')
set(gca, 'FontSize',10);
legend('\sigma_d =0.25', '\sigma_d= 0.5', '\sigma_d =0.75', '\sigma_d = 1')

[David Harra]

Perfect. Thank you so much

Answer 2

[Dyuman Joshi]

You will have to explicitly change the function handle to change its definition

P=pi;
fun = @(x) P*x;
P=3;
%fun(3) is not equal to 3*3=9
fun(3)

ans = 9.4248

L=0:100;
L_bar=25;
sigma_d =[0.25, 0.5, 0.75, 1.0]; 
mu = log(L_bar) - 0.5*(sigma_d).^2;
L_tilda = exp(mu);

% eta_r = P.*exp(-r./L);
%  I want to integrate eta from zero to infinity where P has changing
%  variables to give 4 different plots.


r = 0:100;

% I was just testing below to see if a loop would run before attempting to
% plot anything

for ii=1:numel(sigma_d)
P= 1./(L.*sqrt(2*pi)*sigma_d(ii)).*exp(-log(L/L_tilda(ii)).^2/(2*sigma_d(ii)^2));
fun = @(x) P.*exp(-x./L);
eta_r=integral(fun,0,Inf,'ArrayValued',true)
end

Warning: Inf or NaN value encountered.

eta_r = 1×101

NaN 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0006 0.0030 0.0109 0.0307 0.0718 0.1437 0.2535 0.4022 0.5842 0.7870 0.9943 1.1886 1.3546 1.4810 1.5614 1.5946 1.5834 1.5336 1.4531 1.3500 1.2326

Warning: Inf or NaN value encountered.

eta_r = 1×101

NaN 0.0000 0.0000 0.0003 0.0023 0.0097 0.0269 0.0572 0.1019 0.1598 0.2281 0.3028 0.3800 0.4560 0.5276 0.5924 0.6491 0.6965 0.7345 0.7630 0.7827 0.7940 0.7979 0.7951 0.7867 0.7733 0.7560 0.7354 0.7122 0.6871

Warning: Inf or NaN value encountered.

eta_r = 1×101

NaN 0.0002 0.0060 0.0263 0.0626 0.1109 0.1656 0.2219 0.2764 0.3268 0.3717 0.4106 0.4433 0.4701 0.4914 0.5076 0.5192 0.5268 0.5309 0.5319 0.5303 0.5265 0.5209 0.5137 0.5053 0.4958 0.4855 0.4746 0.4632 0.4514

Warning: Inf or NaN value encountered.

eta_r = 1×101

NaN 0.0099 0.0513 0.1074 0.1642 0.2156 0.2596 0.2959 0.3252 0.3482 0.3658 0.3789 0.3882 0.3942 0.3977 0.3989 0.3984 0.3963 0.3931 0.3889 0.3839 0.3783 0.3722 0.3658 0.3590 0.3521 0.3450 0.3378 0.3305 0.3233

Since you are dividing by 0 (first element of L), you get a NaN and thus the warning from the integral() solver.

[David Harra]

I am not sure I am following everything you are saying. Your code runs but the plots look nothing like what I am expecting. I have tried to read into the function handles but I find everything very vague. I need to be integrating with respect to L I am not sure in your code why you are using x in the integraion and not r.

I am not sure how to amend your code to match the results I require.

[Dyuman Joshi]

"I need to be integrating with respect to L"

You should have mentioned that before. By the definition of the function handle, I took it as you are integrating w.r.t r

Now as you want to integrate w.r.t L, is there a need to define L=0:100?

Or L is a variable in the definition of P as well?

[David Harra]

Aplogies

So my definition of P is

P= 1./(L.*sqrt(2*pi)*sigma_d).*exp(-log(L/L_tilda).^2/(2*sigma_d^2))

So L is definied within P and L defined as a length. I think I could maybe make this a fixed constant of L=100 and keep r=0:100

The integral I am trying to calculate is

I am changing values of L_tilda and sigma 4 times to get 4 different plots.

Hopefully this is more clear :)