Problem with function_handle

I have this code:

rho = 0.04; 
c = 25;
eta = 0.02;
mu = 40;
sigma  = 0.5;
alpha = -rho/eta;

F = @(t) t.^(-alpha-1).*exp(-t.^2./2+((p-mu)./sigma.*sqrt(2.*eta)).*t); 
F_d = @(t)t.^(-alpha-1).*(t.*sqrt(2.*eta))./sigma.*exp(-t.^2./2+((p-mu)./sigma.*sqrt(2.*eta)).*t); 
D = @(p)(exp(((p-mu)./sigma.*sqrt(2.*eta)).^2./4)./gamma(-alpha)).*quadgk(F,0,inf); 
D_d = @(p)(exp((((p-mu)./sigma.*sqrt(2.*eta)).^2)./4)/gamma(-alpha)).*quadgk(F_d,0,inf)+((eta.*(p-mu))./(4.*(sigma).^2)).*(exp((((p-mu)./sigma.*sqrt(2.*eta)).^2)./4)./gamma(-alpha)).*quadgk(F,0,inf); %Funzione cilindrica derivata
psi = @(p)exp((eta.*(p-mu).^2)./(2.*(sigma).^2)).*D; 
psi_d = @(p)(eta.*(p-mu))./(sigma.^2).*exp((eta.*(p-mu).^2)./(2.*(sigma).^2)).*D + exp((eta.*(p-mu).^2)./(2.*(sigma).^2)).*D_d;
Theta = @(p)(p-c).*psi_d-psi;
%
FreeBoundary = fzero(Theta,10)

but it shows me this error

Error using fzero (line 289)
FZERO cannot continue because user supplied function_handle ==> @(p)(p-c).*psi_d-psi
 failed with the error below.

   Undefined function 'times' for input arguments of type 'function_handle'.

Error in prova_2 (line 24)
FreeBoundary = fzero(Theta,10)

Can you tell me why? Where is the error?

Answer 1

[Steven Lord]

You can't add/subtract/multiply/divide function handles so this would not work.

 s = @sin;
 c = @cos;
 sc = @(x) s.*c;
 thisWillError = sc(0:0.1:2*pi);

What you can do is add/subtract/multiply/divide the values returned by function handles instead.

 s = @sin;
 c = @cos;
 sc = @(x) s(x).*c(x);
 thisWillWork = sc(0:0.1:2*pi);

Answer 2

[Stephen23]

You define D as a function:

D = @(p)(exp(((p...

but then later you use it as if it were a numeric (number):

psi = @(p)exp(...).*D;

A function handle is not a numeric, it cannot be multiplied. To fix this you could evaluate that function, which returns a numeric value and this can be multiplied:

psi = @(p)exp(...).*D(p);

[FrancescoMi]

But I need to maintain p as a variable and solve the final equation Theta. I've tried in this way:

rho = 0.04; 
c = 25;
eta = 0.02;
mu = 40;
sigma  = 0.5;
alpha = -rho/eta;

syms p t 
F = t^(-alpha-1)*exp(-t^2/2+((p-mu)/sigma*sqrt(2*eta))*t); 
F_d = t^(-alpha-1)*(t*sqrt(2*eta))/sigma*exp(-t^2/2+((p-mu)/sigma*sqrt(2*eta))*t); 
D = (exp(((p-mu)/sigma*sqrt(2*eta))^2/4)/gamma(-alpha))*int(F,t,0,inf); 
D_d = (exp((((p-mu)/sigma*sqrt(2*eta))^2)/4)/gamma(-alpha))*int(F_d,t,0,inf)+((eta*(p-mu))/(4*(sigma)^2))*(exp((((p-mu)/sigma*sqrt(2*eta))^2)/4)/gamma(-alpha))*int(F,t,0,inf);
psi = exp((eta*(p-mu)^2)/(2*(sigma)^2))*D; 
psi_d = (eta*(p-mu))/(sigma^2)*exp((eta*(p-mu)^2)/(2*(sigma)^2))*D + exp((eta*(p-mu)^2)/(2*(sigma)^2))*D_d; 
Theta = (p-c)*psi_d-psi;

FreeBoundary = solve(Theta,p)

But the values of the integral are really big and the integral doesn't work for more values of rho,eta and sigma.

[Stephen23]

My answer works just fine, when you actually try it:

F = @(p,t) t.^(-alpha-1).*exp(-t.^2./2+((p-mu)./sigma.*sqrt(2.*eta)).*t);
F_d = @(p,t) t.^(-alpha-1).*(t.*sqrt(2.*eta))./sigma.*exp(-t.^2./2+((p-mu)./sigma.*sqrt(2.*eta)).*t);
D = @(p)(exp(((p-mu)./sigma.*sqrt(2.*eta)).^2./4)./gamma(-alpha)).*quadgk(@(t)F(p,t),0,inf);
D_d = @(p)(exp((((p-mu)./sigma.*sqrt(2.*eta)).^2)./4)/gamma(-alpha)).*quadgk(@(t)F_d(p,t),0,inf)+((eta.*(p-mu))./(4.*(sigma).^2)).*(exp((((p-mu)./sigma.*sqrt(2.*eta)).^2)./4)./gamma(-alpha)).*quadgk(@(t)F(p,t),0,inf); %Funzione cilindrica derivata
psi = @(p)exp((eta.*(p-mu).^2)./(2.*(sigma).^2)).*D(p);
psi_d = @(p)(eta.*(p-mu))./(sigma.^2).*exp((eta.*(p-mu).^2)./(2.*(sigma).^2)).*D(p) + exp((eta.*(p-mu).^2)./(2.*(sigma).^2)).*D_d(p);
Theta = @(p)(p-c).*psi_d(p)-psi(p);

And tested:

>> Theta(1.2)
ans =   7.4491e+051
>> Theta(-pi)
ans =   1.8295e+064

You might like to plot the data before trying to use fzoer, because it does not seem to have a zero crossing.