Double integral where one limit can only be solved numerically
1 次查看(过去 30 天)
显示 更早的评论
Hello All, This is my first question in this forum.
I'm trying to do a double integral using qua2d function.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/158527/image.png)
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,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/158528/image.png)
Any help setting up this integral would be much appreciated.
Thank You.
3 个评论
Star Strider
2014-3-15
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)?
采纳的回答
Roger Stafford
2014-3-14
编辑:Roger Stafford
2014-3-14
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);
更多回答(1 个)
Roger Stafford
2014-3-14
编辑:Roger Stafford
2014-3-14
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!