3

I've got a question about the implementation of a double integral in MATLAB.

It's known that enter image description here

Making use of

k1 = 1E-04:0.001:1E+04; k2 = 1E-04:0.001:1E+04; k3 = 1E-04:0.001:1E+04;

the above procedure(calling the formula of F11,F22 and F33) leads to the results shown below: enter image description here

Now comes the question:

I would like to achieve the same results by means of a double integral (or nested single integral in k2 and k3) involving only phi11,phi22 and phi33, hence without calling directly the formula for F11, F22 and F33.

Up to now I'm using this code (restricted to F11 calculation):

clc,clear all,close all
k1 = (1E-04:0.01:10000);
k2 = (1E-04:0.01:10000);
k3 = (1E-04:0.01:10000);
F_11_benchmark = ((9/55)*1.453)*(1./(1+k1.^2).^(5/6));
count = 0;
F_11 = zeros(1,numel(k1));
for i = 2:numel(k1)
count = count + 1;
phi_11 = @(k2,k3) (1.453./(4.*pi)).*((k2.^2+k3.^2)./((1 + k1(count).^2 + k2.^2+k3.^2).^(17/6)));
F_11(count) = dblquad(phi_11,k2(i-1),k2(i),k3(i-1),k3(i));
end

not leading to the same benchmark results shown in the last figure.

How would you suggest to tackle this problem? I thank you in advance.

Best regards, FPE

Community
  • 1
  • 1
fpe
  • 2,700
  • 2
  • 23
  • 47

1 Answers1

1

You are using dblquad wrong. The way you are doing it calculates for each k1 the integral over the little box [k1(i-1) k1(i)]x[k1(i-1) k1(i)], which in no way approximates the whole domain of integration. Use the integral2 command instead: try just using

F_11(count) = integral2(phi_11,-100,100,-100,100);

or

F_11(count) = integral2(phi_11,-Inf,Inf,-Inf,Inf);

which should be more accurate but will be slower. You can probably also get away with using a much larger resolution for k1. I used

k1=10.^(-4:.01:3);

to sort of match your graph.

RussH
  • 329
  • 2
  • 17
  • It works, but I'm not sure if this procedure is the one mathematically rigorous. Within the scientific paper I've been reading looks more one has to integrate on small control volumes. Btw, thanks for the hint, was precious. – fpe Jan 02 '13 at 15:51
  • could you post the methods you believe would let the procedure being more rigorous? thanks in advance. – fpe Jan 02 '13 at 19:30
  • @fpe dblquad takes care of that for you. It iterates on different refinements of the box you define. the command i mentioned evaluates the integral to a specified tolerance. it takes care of all the refinement for you. you can make it more accurate by increasing the domain of integration from [-100,100]x[-100,100] to something bigger . Additionally you can change the tolerance dblquad uses (look up the help file for that function). Also look at integral2, you can use infinite limits with that one. Don't worry about the rigor because MATLAB has already worried about it for you :) – RussH Jan 02 '13 at 21:28
  • Actually I would strongly recommend using integral2 instead of dblquad. Same arguments, editing answer – RussH Jan 02 '13 at 21:29
  • I had a look at the dblquad help and read all the features this function embodies. A 1E-06 tolerance is already fine enough, so that I'm just gonna use the procedure you advise. – fpe Jan 02 '13 at 21:30
  • you're again right: integral2 is retrieving better and more accurate results than dblquad. – fpe Jan 02 '13 at 22:00