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.

[Star Strider]

What are typical values for p, q, r, and V? Can lambda be complex?

[Sam]

Thank you for your response Star Strider. lambda cannot be complex.p and q can be any number. let's say -1<p,q<5 and 0<r<1. 0<V<5

[Star Strider]

I’ve been following this discussion out of interest.

Consider writing this up when you’re finished with it as a function file or series of function files and submitting it to the File Exchange.

I’d also like to experiment with it on my own when I have the time. What function are you using for f(x,y) (the denominator of your integrand)? If that’s a topic of your research you’d prefer not to reveal at present, what is a typical f(x,y)?

Answer 1

[Roger Stafford]

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

[Sam]

Thank You Roger,

No, I do not anticipate any singularities by the division of f(x,y).

I didn't know I could use a function call as one of the limits in the 'integral2'. I used your code and edited it to my needs. However it gave me an answer, 6 orders of magnitudes off. I'm suppose to get values in the range of 10^-3, but it returned a value in the range 10^3. I couldn't find any mistake in the edited code. However now I have an idea on how would I go about doing something like this. I appreciate your help and hopefully I'll be able to find and fix the problem. I have just one question, I couldn't find any relationship between the limits of the outer integral and 'X'(input of the function). Could you please explain me what values are passed to the 'samymax' function when you call it inside 'integral2'. Thank you.

[Roger Stafford]

I meant to suggest this to you before. Before trying to do the double integration you should thoroughly test this method of using 'fzero' to find the value of 'lambda' properly. In place of the call on 'integral2' you should place some routine in the same code I gave you that calls on 'samymax' with a number of reasonably realistic vectors of x values to see if everything is working all right with it. This routine, for each x and the corresponding lambda it receives from 'samymax', would compute the quantity

   p-q-lambda-sqrt(r*lambda+exp(lambda-x)

to see if its value is consistently zero (within round-off error limits of course.) You should also test that the parameters p, q, r, and V are being passed correctly to 'samymax' by displaying them from within.

In other words I am urging you not to give up on this effort too easily. Sometimes there is only some trivial mistake in such code.

In answer to your question "Could you please explain me what values are passed to the 'samymax' function", only Mathworks knows the intimate details of their calls to obtain the variable integration limits on the inner integral, but presumably it is a vector of x values covering some range of x in the outer integration process. It is apparently preparing to cover some small range of x with a number of inner integrations, and it needs the limits for each one all at once before it can undertake that. For accuracy purposes it may have to repeat this over the same range of x for a finer resolution of x's. I don't think you need to concern yourself with how it's done, however, Your job is simply to pass to it the correct upper limits for a batch of x's it passes to you. That is why I think you should concentrate on making sure 'samymax' is doing its thing properly.

I wrote the code to allow for a matrix of values being passed to 'samymax' because of Mathworks' statement "The functions ymin=c(X) and ymax=d(X) must accept matrices and return matrices of the same size with corresponding values" made in the documentation for 'quad2d'. However, I seriously doubt that they would actually send a two-dimensional matrix for inner integration limits, but there is no particular penalty in providing for that just in case it turned out to be true.

[Sam]

X is always a 14 element vector consisting values between upper and lower limits of the outer integral. It doesn't pass elements starting from the lower limit to the upper limit. For instance, when I set the upper limit to 3, the last set of elements passed to X were, "0.6974, 0.6971, 0.6964, 0.6956, 0.6947, 0.6941, 0.6937, 0.6935, 0.6932, 0.6925, 0.6917, 0.6908, 0.6902, 0.6898". Before that a set of values in the range of "2.7xx" were passed. I didn't get a time to work on this yesterday. I will work on it today and I'll be sure to let you know how it goes. Thank you Roger.

[Sam]

Turns out I had a small typo in my integral. Now it returns the correct answer. Thank You so much Roger.

Answer 2

[Roger Stafford]

It's good that you corrected that error replacing V with x. That's an important distinction in the integration process. However, it doesn't really change the difficult part of your problem, and that is the solution to your equation

 p-q-lambda = sqrt(r*lambda+exp(lambda-x))

That is where you need to concentrate your energies.

As for the variable upper limit of integration, namely, a lambda which varies with x, both 'quad2d' and 'integral2' provide for the limits of integration on the inner integral to be variable in accordance with function handles. In other words they are able to integrate over a non-rectangular x,y region. You should read up the description of these function handles carefully to abide by their rules. That fact that your upper limit is not in the form of an explicit expression does not matter. If you can write a function that has the ability to return lambda's that correspond to entered x's, that is all that is needed to do the integration with either of these integration functions.

If I were doing your problem I would be tempted to use matlab's 'fzero' function to produce 'lambda'. Obviously you cannot do it with an anonymous function. It has to be an m-file function or subfunction. The integration routine is going to call on it with x values in the form of either vectors or perhaps matrices. You have to write the function to deal with each x value as a separate problem probably using for-loops and produce a corresponding lambda and then return with a like vector or matrix of these lambda values.

Within these for-loops you are faced with the same problem as before - given a single value of x, find the lambda value that satisfies the above equation. The 'fzero' routine allows you to enter a two-valued vector as your "estimate" argument, which I recommend. As the lower value you should use the least possible value that lambda could have. It should be so low that it always satisfies the inequality

p-q-lambda - sqrt(r*lambda+exp(lambda-x)) > 0

Perhaps zero would always serve there but you should check that out. It should also satisfy

r*lambda+exp(lambda-x) >= 0

so as to avoid complex results which would be disturbing to 'fzero'.

Choosing the upper value of the interval is perhaps a little trickier. The upper value must reverse the inequality above so that

p-q-lambda - sqrt(r*lambda+exp(lambda-x)) < 0

To accomplish this you may have to do a little iteration of your own, say repeatedly doubling some starting value until the inequality holds. Or perhaps you might discover some more sophisticated method. Just make sure it doesn't take any longer than the iteration done by 'fzero' does.

Once you have these two inequalities both satisfied, your solution will be pinned somewhere between them and 'fzero' will certainly "zero" in on the solution, and rather rapidly if I am not mistaken.

You will have to face the fact that computation will require more time than would be usual for a double integration procedure because of the heavy demands made on 'fsolve' to find 'lambda'.

That is an outline of how I would do the problem. I hope you manage to make it work.

[Sam]

Hello Roger, this is what I've been looking for. After reading this, I have a rough idea what I'm suppose to do. but I'm a little confused. I'm not sure how would I use an m-file or a sub-function to call the integration and how would I use for loops to solve for lambda at each x in the integration. can you show me using a simple example or can you show me using simple code (even pseudo code) on what exactly would you have done to solve this. Thank You.