How can I use a numerical solution to an equation as an anonymous function for an integration?

I am attempting to do a numerical triple integral using anonymous functions to set up the integrand. Up until now, the approximations I have made in deriving the equations have been sufficient, and to do the integral I have the code:

kf = @(k0,psi) k00+k06*cos(6*(psi))+k012*cos(12*(psi))+k31*sin(c*k0/N).*cos(3*(psi))+k20*cos(2*c*k0/N);

vperp = @(k0,psi) -3*hbar/m.*(2*k06.*sin(6.*(psi))+4*k012.*sin(12.*(psi))+k31.*sin(c.*k0./N).*cos(3.*(psi))); 

v = @(k0,psi) hbar*kf(k0,psi)/m;

vx = @(k0,psi) cos(psi).*v(k0,psi)-sin(psi).*vperp(k0,psi); 

integrand = @(k0,psi,psipr) vx(k0-k00*cos(phi)*sin(theta),psi).*vx(k0-k00*cos(phi)*sin(theta),psipr).*exp(psipr.*cos(theta)./wt0);

limit = @(k0,psi) psi*cos(theta)+phi;

f = e^2*cos(theta)*m/(4*pi^3*hbar^2*omega_c)*integral3(integrand,-pi/c,pi/c,phi,2*pi/cos(theta)+phi,-Inf,limit);

phi is an angle between 0 and 2pi, N=3, and hbar and m are Planck's reduced constant and the mass of the electron respectively.

However, now I have changed some parameters so that my original approximations are incorrect, so now the 1st argument of integrand must be, instead of that shown above the solution to the equation:

*x* = k0 - (k00 + k31*sin( *x* *c/3)*cos(3*psi))*cos(phi)*sin(theta

This can only be solved numerically, obviously.

I would love to know how to get this to work, either using the method I've currently been using, or by another method, so if anyone has any ideas, that would be great!

Thanks

Answer 1

[Matt J]

 integrand=@(k0,psi,psipr)...
           fzero(@(x)k0-(k00+k31*sin(c/3)*cos(3*psi))*cos(phi)*sin(theta)-x,...
            [xlow,xhigh] )

[Joseph]

Thank you very much, this looks very useful! I presume that xlow and xhigh are the endpoints of the interval in which the solution is to be found? There is one further problem I'm having though - I expect the solution to be near k0-k00*cos(phi)*sin(theta), but if I put the end points of the interval at k0-k00*cos(phi)*sin(theta) +k31 and -k31, I get an error about the second argument not being a scalar or vector of length 2. I'm not really sure how to fix this, which seems to come from the fact that the endpoints are also functions of k0. Do you have any further ideas? I would really appreciate it. Thanks for your help so far!

[Matt J]

Ah, I see. You will have to write a function that accepts arrays k0,psi,psipr and solves the equation for each combination of input values in a loop.

[Mike Hosea]

If you have a function f(x,y,z) that works correctly when x, y, and z are scalars, then you can use

fv = @(x,y,z)arrayfun(f,x,y,z)

to handle the looping over each element when x, y, and z are non-scalar.

[Joseph]

Thank you for your help so far. I've tried putting in arrayfun, but it seems to be having problems with it. I'm pretty certain I have a function that solves the equation given some scalar values for the variables, but I'm getting an error corresponding to this line:

sigxxintegrand = @(k0,psi,psipr) arrayfun(integrand,k0,psi,psipr);

integrand is a function that has 3 input variables, which should be k0, psi and psipr. However, when I try running the code, I get the error 'Not enough input arguments' in this line. How am I misusing arrayfun?

My full code is:

limit = @(psi) psi.*cos(theta)+phi;

function result = integrand(k0,psi,psipr)
  kf = @(x) k00+k06*cos(6*(psi))+k012*cos(12*(psi))+k31*sin(c.*x./N).*cos(3*(psi))+k20*cos(2*c.*x./N);
  kfpr = @(x) k00+k06*cos(6*(psipr))+k012*cos(12*(psipr))+k31*sin(c*x/N).*cos(3*(psipr))+k20*cos(2*c*x/N); 
  vperp = @(x) -3*hbar/m.*(2*k06.*sin(6.*(psi))+4*k012.*sin(12.*(psi))+k31.*sin(c.*x./N).*cos(3.*(psi))); 
      vperppr = @(x) -3*hbar/m.*(2*k06.*sin(6.*(psipr))+4*k012.*sin(12.*(psipr))+k31.*sin(c.*x./N).*cos(3.*(psipr))); 
      v = @(x) hbar*kf(x)/m;        
      vpr = @(x) hbar*kfpr(x)/m;        
      vx = @(x) cos(psi).*v(x)-sin(psi).*vperp(x); 
      vxpr = @(x) cos(psipr).*vpr(x)-sin(psipr).*vperppr(x); 
      fzerolowerlimit = k0-k00*cos(phi)*sin(theta)-k31; %Limits for fzero function
      fzeroupperlimit = k0-k00*cos(phi)*sin(theta)+k31;
      limits = [fzerolowerlimit,fzeroupperlimit];
      equation = @(kz) k0-kz-k00*cos(phi)*sin(theta)-k31*sin(kz*c/3)*cos(3*psi)*cos(phi)*sin(theta); %Equation to be solved by fzero function
      argument = fzero(@(kz)equation(kz),limits); %Defines the equation that must be solved for kz and then solves it, giving the output as the argument required
  result = vx(argument).*vxpr(argument).*exp(psipr.*cos(theta)./wt0); %Defines the integrand
end

sigxxintegrand = @(k0,psi,psipr) arrayfun(integrand,k0,psi,psipr);

Sxx1 = @(k0,psi) integral(@(psipr)sigxxintegrand(k0,psi,psipr),-Inf,limit(psi),'ArrayValued',true);
Sxx2 = @(k0) integral(@(psi)Sxx1(k0,psi),phi,2*pi/cos(theta)+phi,'ArrayValued',true);
Sxx = e^2*cos(theta)*m/(4*pi^3*hbar^2*omega_c)*integral(@(k0)Sxx2(k0),-pi/c,pi/c,'ArrayValued',true);

[Matt J]

Can't run your full code since data for e, theta, m, etc... is not provided. However, I do not get an error at the line you've shown. It does make it to the final line

 Sxx = e^2*cos(theta)*m/(4*pi^3*hbar^2*omega_c)*integral(@(k0)Sxx2(k0),-pi/c,pi/c,'ArrayValued',true);

[Joseph]

Ah, apologies, the values I am using are:

phi=10*pi/180;
theta=-50*pi/180;
k00=1e10;
k06=0.04e10;
k012=0.007e10;
k20=-0.0025e10;
k31=0.1e10;
hbar=1.05e-34;
m=9.109e-31;
N=3;
c=17.837e-10;

That does seem a bit odd though, I've tried running it again, exactly as I wrote (although I removed the list of constants in front of the integral in the last line, as they aren't really necessary), with those values, and I still get the error of there not being enough input values. Could it be something to do with the version of Matlab? Again, thank you for your help.

[Matt J]

You're missing an @ operator. You need to specify "integrand" using a function handle, not just its name.

 sigxxintegrand = @(k0,psi,psipr) arrayfun(@integrand,k0,psi,psipr);

[Joseph]

Ahh thank you, I thought it would end up being something stupid like that. I am still experiencing problems with the code, as it gives me NaN for the final result, which seems to be an issue with the integration, as the arrayfun part seems to be working fine. I believe the issue is with the limit which is a function of psi, as if I replace it with 0 the program seems to work fine. I have also slightly changed the definition of limit so that we have

function l = limit(k0,psi)
    l = psi.*cos(theta)+phi;
end

Then we have

integrallimit = @(k0,psi) arrayfun(@limit,k0,psi);
Sxx1 = @(k0,psi) integral(@(psipr)sigxxintegrand(k0,psi,psipr),-Inf,integrallimit(k0,psi),'ArrayValued',true);

If no immediately obvious reason or solution to this springs to mind, please don't worry, you've been a great help, thank you very much!