MATLAB: Double integral where one limit can only be solved numerically

Hello All, This is my first question in this forum.
I'm trying to do a double integral using qua2d function.
Here, the upper limit of the inner integral lambda is a function of 'V' which is the variable of the outer integral. But the problem is that there's no explicit equation for lambda. lambda can only be solved numerically using an equation of the form,
Any help setting up this integral would be much appreciated.
Thank You.

Best Answer

The following is approximately what I think you need to do. It may need some adjusting here and there. I am unable to test it directly on my own primitive version of matlab, but at least it will give you an idea of what I was trying to say in the outline. Note that I have assumed that the value zero would be all right as a lower limit in the "estimate" interval passed to 'fzero'. It should be okay if p-q>=1. Otherwise setting it properly would be a more difficult task. The upper end of the interval, 'up', is doubled until the inequality reverses. I also assume r>=0 to ensure that p-q-lambda-sqrt(r*lambda+exp(lambda-x) decreases as lambda increases. I am assuming the parameters p,q, r, and V can be passed to the subfunction 'samymax' in the manner I've shown here. Let me know if you encounter unsurmountable difficulties with this. Note that I haven't placed any options as to error tolerances in either 'fzero' or 'integral2'. You may want to specify these to guarantee good accuracy.
I do have one question. Do you anticipate encountering any singularities in the integrand from the division by f(x,y) for the specified ranges of x and y? Matlab's 'integral2' function is supposed to be able to handle singularities in the integrand but only if they are of an integrable nature- that is, only if the integral isn't divergent.
Here is my suggested code:
 % The parameters
 p = ... % We assume p-q >= 1
 q = ...
 r = ...
 V = ...
 % The integrand as an anonymous function of x and y
 samint = @(x,y)A*exp(y-x)./f(x,y); % <-- Fix this up so it's correct
 % The inner integral upper limit as a function of x
 function lambda = samymax(X)
   [m,n] = size(X);
   lambda = zeros(m,n);
   for ix1 = 1:m
     for ix2 = 1:n
       x = X(ix1;ix2);
       up = 1; % We assume r >= 0
       while p-q-up-sqrt(r*up+exp(up-x) >= 0
         up = 2*up;  % Double 'up' each time until inequality is "<"
       end
       pqrx = @(lam)p-q-lam-sqrt(r*lam+exp(lam-x);
       % Assume p-q>=1, so that lower interval value of zero is valid 
       lambda(ix1,ix2) = fzero(pqrx,[0,up]);
     end
   end
 end
 % The double integral itself
 I = integral2(samint,0,V,0,@samymax);
Related Question