Double Integration problem

I am trying to integrate the following function.This is my code in matlab and the error messages:

pc=0;

f1 = (log(x) - log(31.1)./0.48).^2;

f2 = 2 .* pc .* ((log(x) - log(31.1)) ./ 0.48) .* ((log(y) - log(4.5)) ./ 0.58);

f3 = ((log(y) - log(4.5)) ./ 0.58) .^ 2;

P = @(x,y)(exp((-0.5 .* (f1 - f2 + f3)) ./ (1 - pc .^ 2))). / (2 .* pi .* x .* y .* 1.21 .* 0.9 .* sqrt(1 - pc .^ 2));

Q = dblquad(P,10,10000,10,10000)

 Error using ==> mldivide
Matrix dimensions must agree.

Error in ==> @(x,y)(exp((-0.5.*(f1-f2+f3))./(1-pc.^2)))/(2.*pi.*x.*y.*1.21.*0.9.*sqrt(1-pc.^2))

Error in ==> quad at 76 y = f(x, varargin{:});

Error in ==> dblquad>innerintegral at 77 Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});

Error in ==> quad at 76 y = f(x, varargin{:});

Error in ==> dblquad at 53 Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

I probably have syntax errors or maybe a dot is misplaced.I am new in Matlab and any help would be highly appreciated.Thank you.

Answer 1

[the cyclist]

Does this do what you expect? I took care of the syntax error that Walter pointed out, as well as the function specification I mentioned in my other answer now.

pc=0;

f1 = @(x)((log(x) - log(31.1)./0.48).^2);

f2 = @(x,y)(2 .* pc .* ((log(x) - log(31.1)) ./ 0.48) .* ((log(y) - log(4.5)) ./ 0.58));

f3 = @(y)(((log(y) - log(4.5)) ./ 0.58) .^ 2);

P = @(x,y)((exp((-0.5 .* (f1(x) - f2(x,y) + f3(y))) ./ (1 - pc .^ 2)))./ (2 .* pi .* x .* y .* 1.21 .* 0.9 .* sqrt(1 - pc .^ 2)));

Q = dblquad(P,10,10000,10,10000)

[Geo]

Thanks a lot for the answer. You are my "savior". I would like to ask something else . Can i use the dblquad function with infinity or -infinity?? I mean something like this:

Q = dblquad(P,10,Inf,10,Inf) or

Q = dblquad(P,-Inf,10,10,Inf). Or is there another way to approximate the value of infinity??

[the cyclist]

I don't think MATLAB is going to be able to figure out a numerical integral with an Inf limit like that. I can think of a few possible approaches. (1) Try a change of variable, such that the limits become finite; (2) Try to infer the limit of the integral by trying larger and larger intervals [not rigorous, and probably not going to work]; (3) Try using the Symblic Math Toolbox.

[Also if you found this answer helpful, it is best to "accept" it, for future users to see it might be most helpful]

[Walter Roberson]

http://www.mathworks.com/help/techdoc/ref/quadgk.html supports infinite limits, as does http://www.mathworks.com/help/techdoc/ref/quad.html

See also the discussion http://www.mathworks.com/matlabcentral/newsreader/view_thread/134595

and it would probably be good to read the discussion of methods in

http://www.mathworks.com/help/techdoc/ref/quadl.html

Note that dblquad() allows you to specify which method to call, so you *might* be able to specify quadgk as the method with infinite limits, but that is not something I have ever tried.

[Geo]

Thanks again for all your help and the time spent.Ι probably use as an upper limit 10^16 or something like that.

Answer 2

[Walter Roberson]

In your statement

P = @(x,y)(exp((-0.5 .* (f1 - f2 + f3)) ./ (1 - pc .^ 2))). / (2 .* pi .* x .* y .* 1.21 .* 0.9 .* sqrt(1 - pc .^ 2));

notice that you have ). / when instead you want ) ./

Answer 3

[the cyclist]

I don't get that error, because I cannot get past the syntax error in the definition of P. That line, as you have placed the parentheses, gives me an error: "Unexpected MATLAB operator."

Also, I assume f1,f2, and f3 are supposed to be their own function definition, but they are not, as you have written them.