Double integral where one limit can only be solved numerically

3 views (last 30 days)
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.
  3 Comments
Sam
Sam on 13 Mar 2014
Edited: Sam on 13 Mar 2014
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
Star Strider on 15 Mar 2014
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)?

Sign in to comment.

Accepted Answer

Roger Stafford
Roger Stafford on 14 Mar 2014
Edited: Roger Stafford on 14 Mar 2014
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);
  4 Comments
Sam
Sam on 16 Mar 2014
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
Sam on 16 Mar 2014
Turns out I had a small typo in my integral. Now it returns the correct answer. Thank You so much Roger.

Sign in to comment.

More Answers (1)

Roger Stafford
Roger Stafford on 14 Mar 2014
Edited: Roger Stafford on 14 Mar 2014
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.
  1 Comment
Sam
Sam on 14 Mar 2014
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.

Sign in to comment.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Tags

No tags entered yet.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!