Search code examples
matlabintegral

Benchmark integral vs. Double integral procedure

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

Solution

  • 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.