numerical integration 2d function of three variables and plotting

Hey guys.

The problem is the following: I have a (quiet complicated) function f(x,y,z) and want to integrate it over x and y to later plot the result versus z. I have seen the existing question about getting a symbolic function of z after the two integrations but that does not help me: my function is not separable. I also don't need an explicit function of z as I just want to plot the 2d integral against z.

I tried "int" which gives the error: undefined function 'int' for input arguments of type 'function handle';

"quad2d" giving: A must be a finite, scalar, floating point constant; and "integral2" giving me millions of errors.

Any ideas? I would really appreciate it!

clf reset;
Delta_inf=1.227;
omega = 0:1/100:7*Delta_inf;
Delta_i = 1.35;
beta = 2.00;
epsilon_f = 9479;
f_R1 = @(epsilon_1) 1./((Delta_inf-epsilon_1).*sqrt((epsilon_1-epsilon_f).^2+Delta_i.^2));
f_R2 = @(epsilon_2) 1./((-Delta_inf-epsilon_2).*sqrt((epsilon_2-epsilon_f).^2+Delta_i.^2));
f_1 = integral(f_R1,-Inf,Inf);
f_2 = integral(f_R2,-Inf,Inf);
syms epsilon_1 epsilon_2
E_1 = sqrt(epsilon_1.^2+Delta_inf.^2);
E_2 = sqrt(epsilon_2.^2+Delta_inf.^2);
N_1 = ((epsilon_1-epsilon_f)*f_1+beta).^2+(Delta_i*f_1).^2+pi.*pi;
N_2 = ((epsilon_2-epsilon_f)*f_2+beta).^2+(Delta_i*f_2).^2+pi.*pi;
gamma_1 = -sqrt(1-(N_1-sqrt(N_1.^2-(2.*pi.*Delta_i.*beta./E_1).^2))./(2.*pi.*pi));
gamma_2 = -sqrt(1-(N_2-sqrt(N_2.^2-(2.*pi.*Delta_i.*beta./E_2).^2))./(2.*pi.*pi));
resa = @(epsilon_1, epsilon_2) (gamma_1+gamma_2).*(1-Delta_inf./(E_1.*E_2)).*(dirac(omega+E_1+E_2)-dirac(omega-E_1-E_2));
intresa = quad2d(resa,-Inf,Inf,-Inf,Inf);
resb = @(epsilon_1,epsilon_2) (gamma_1-gamma_2).*(1+Delta_inf./(E_1.*E_2)).*(dirac(omega+E_1-E_2)-dirac(omega-E_1+E_2));
intresb = quad2d(resb,-Inf,Inf,-Inf,Inf);
Realsigmaa = 1./(8.*omega).*intresa;
Realsigmab = 1./(8.*omega).*intresb;
xaxis = omega./Delta_inf;
Realsigma = Realsigmaa+Realsigmab;
plot(xaxis, Realsigma);

the quad2d error:

Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
Error in untitled (line 27)
intresa = integral2(resa,-Inf,Inf,-Inf,Inf);

The int error:

Undefined function 'int' for input arguments of type 'function_handle'.
Error in Versuch (line 28)
intresa = int(resa,-Inf,Inf,-Inf,Inf);

the integral2 error:

Error using integralCalc/finalInputChecks (line 522)
Input function must return 'double' or 'single' values. Found 'sym'.
Error in integralCalc/iterateScalarValued (line 315)
                finalInputChecks(x,fx);
Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);
Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...
Error in
integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x))
(line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...
Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);
Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);
Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);
Error in integral2Calc (line 7)
    [q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);
Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
Error in Versuch (line 28)
intresa = integral2(resa,-Inf,Inf,-Inf,Inf);

[Jan]

Please post the code and copies of the error messages. It is hard to suggest a fix for "millions of errors".

[Walter Roberson]

quad2d cannot be used for infinite ranges.

Answer 1

[Walter Roberson]

The integral of f_R1 over -inf to +inf is undefined because of the singularity at +1.227 . You might as well replace it with a random number because you will get some nonsense value when you integral() it.

The integral of f_R2 over -inf to +inf is undefined because of the singularity at -1.227 . You might as well replace it with a random number because you will get some nonsense value when you integral() it.

dirac() is strictly symbolic and cannot be used with quad2d or integral2()

omega is a vector, so your resa calculates a vector. vectors cannot be integrated with integral2

Delta_inf=1.227;
omega = 0:1/100:7*Delta_inf;
Delta_i = 1.35;
beta = 2.00;
epsilon_f = 9479;
f_R1 = @(epsilon_1) 1./((Delta_inf-epsilon_1).*sqrt((epsilon_1-epsilon_f).^2+Delta_i.^2));
f_R2 = @(epsilon_2) 1./((-Delta_inf-epsilon_2).*sqrt((epsilon_2-epsilon_f).^2+Delta_i.^2));
%f_1 = integral(f_R1,-Inf,Inf);
%f_2 = integral(f_R2,-Inf,Inf);
disp('randomizing f_1 since it is undefined')
f_1 = randn()
disp('randomizing f_2 since it is undefined')
f_2 = randn()
syms epsilon_1 epsilon_2
E_1 = sqrt(epsilon_1.^2+Delta_inf.^2);
E_2 = sqrt(epsilon_2.^2+Delta_inf.^2);
N_1 = ((epsilon_1-epsilon_f)*f_1+beta).^2+(Delta_i*f_1).^2+pi.*pi;
N_2 = ((epsilon_2-epsilon_f)*f_2+beta).^2+(Delta_i*f_2).^2+pi.*pi;
gamma_1 = -sqrt(1-(N_1-sqrt(N_1.^2-(2.*pi.*Delta_i.*beta./E_1).^2))./(2.*pi.*pi));
gamma_2 = -sqrt(1-(N_2-sqrt(N_2.^2-(2.*pi.*Delta_i.*beta./E_2).^2))./(2.*pi.*pi));
resa = @(epsilon_1, epsilon_2) (gamma_1+gamma_2).*(1-Delta_inf./(E_1.*E_2)).*(dirac(omega+E_1+E_2)-dirac(omega-E_1-E_2));
%intresa = quad2d(resa,-Inf,Inf,-Inf,Inf);
intresa = vpaintegral(vpaintegral(resa(epsilon_1,epsilon_2),epsilon_1,-Inf,Inf),epsilon_2,-Inf,Inf)
resb = @(epsilon_1,epsilon_2) (gamma_1-gamma_2).*(1+Delta_inf./(E_1.*E_2)).*(dirac(omega+E_1-E_2)-dirac(omega-E_1+E_2));
%intresb = quad2d(resb,-Inf,Inf,-Inf,Inf);
intresb = vpaintegral(vpaintegral(resb(epsilon_1,epsilon_2),epsilon_1,-Inf,Inf),epsilon_2,-Inf,Inf)
Realsigmaa = 1./(8.*omega).*intresa;
Realsigmab = 1./(8.*omega).*intresb;
xaxis = omega./Delta_inf;
Realsigma = Realsigmaa+Realsigmab;
plot(xaxis, Realsigma);

[Walter Roberson]

the vpaintegral are taking a long time...

You have

Dirac(sqrt(epsilon_1^2+1695074960848869/1125899906842624)+sqrt(epsilon_2^2+1695074960848869/1125899906842624)-1/100)

and

Dirac(sqrt(epsilon_1^2+1695074960848869/1125899906842624)+sqrt(epsilon_2^2+1695074960848869/1125899906842624)+1/100)

Dirac is 0 everywhere except where its input is 0. But for real valued epsilon_1,

sqrt(epsilon_1^2+1695074960848869/1125899906842624)

is a minimum of 3*sqrt(188341662316541)*(1/33554432) which is about 1.227. The sqrt() part after that cannot be negative, so the Dirac() is being asked to operate on something that is at least 1.227 - 1/100 which will always be positive. So the first Dirac() will always output 0.

Likewise, the second Dirac() is being asked to operate on the sum of three non-negative values, and cannot possibly generate 0, so the second Dirac will always be 0 as well.

Therefore the entire difference of dirac term will always be 0, so the integral will come out as 0.

[Walter Roberson]

And after the long vpaintegral, indeed I got out vectors of exact 0 for both functions.

[Rebecca Müller]

sorry maybe I did a huge syntax mistake...I don't want omega = 1/100, but i want omega to take values in the interval [0,7*Delta_inf] and then plot Realsigma against it. So it's true that the first dirac is always 0 for real epsilon, thanks for the advice, but the second one is not always zero, at least if both of the epsilons are smaller than ca. 4.1 (an example would be epsilon_1=4 and epsilon_2=4.1 which would have a dirac peak in omega/Delta_inf~6,9)...

[Walter Roberson]

It doesn't really matter, because your f_R1 and f_R2 are undefined because of the 1/(Delta_inf-epsilon_1) and 1/(-Delta_inf-epsilon_2)

[Rebecca Müller]

Any suggestion how to deal with the singularity? Does matlab know how to deal with cauchy principal values?

[Walter Roberson]

The symbolic toolbox int() has a 'PrincipalValue', true option. Unfortunately it is not strong enough to be able to find the principal values for those two functions. They are

-(2000/89828182862029)*89828182862029^(1/2)*arctanh((9477773/89828182862029)*89828182862029^(1/2))

and

-(2000/89874705794029)*89874705794029^(1/2)*arctanh((9480227/89874705794029)*89874705794029^(1/2))

I revised my code to use those values and to run each of the omega value in a series. In the subset I tested, the values all came out 0, but it is possible that I did not happen to run the series far enough. I am re-running with the full set.

[Walter Roberson]

It took a while, but I confirm that with those principal values, that the integral is 0 for all of the omega values 0:1/100:7*Delta_inf . The exception is for omega 0, which leads to NaN.

Q = @(v) sym(v, 'r');

Delta_inf = Q(1.227);
Delta_i = Q(1.35);
beta = 2.00;
epsilon_f = 9479;
syms epsilon_1 epsilon_2
f_R1 = @(epsilon_1) 1./((Delta_inf-epsilon_1).*sqrt((epsilon_1-epsilon_f).^2+Delta_i.^2));
f_R2 = @(epsilon_2) 1./((-Delta_inf-epsilon_2).*sqrt((epsilon_2-epsilon_f).^2+Delta_i.^2));
f_1 = int(f_R1(epsilon_1), epsilon_1, -Inf, Inf, 'PrincipalValue', true)
f_2 = int(f_R2(epsilon_2), epsilon_2, -Inf, Inf, 'PrincipalValue', true)
magic1 = Q(89828182862029);
f_1 = -Q(2000)*sqrt(magic1)*atanh(Q(9477773)*sqrt(magic1)*(1/magic1))*(1/magic1)
magic2 = Q(89874705794029);
f_2 = -2000*sqrt(magic2)*atanh(Q(9480227)*sqrt(magic2)*(1/magic2))*(1/magic2)

E_1 = sqrt(epsilon_1.^2+Delta_inf.^2);
E_2 = sqrt(epsilon_2.^2+Delta_inf.^2);
N_1 = ((epsilon_1-epsilon_f)*f_1+beta).^2+(Delta_i*f_1).^2+pi.*pi;
N_2 = ((epsilon_2-epsilon_f)*f_2+beta).^2+(Delta_i*f_2).^2+pi.*pi;
gamma_1 = -sqrt(1-(N_1-sqrt(N_1.^2-(2.*pi.*Delta_i.*beta./E_1).^2))./(2.*pi.*pi));
gamma_2 = -sqrt(1-(N_2-sqrt(N_2.^2-(2.*pi.*Delta_i.*beta./E_2).^2))./(2.*pi.*pi));

omega_vals = 0:1/100:7*Delta_inf;
N = length(omega_vals);

xaxis = zeros(1, N, 'sym');
Realsigma = zeros(1, N, 'sym');

for K = 1 : N
    omega = omega_vals(K);
    resa = @(epsilon_1, epsilon_2) (gamma_1+gamma_2).*(1-Delta_inf./(E_1.*E_2)).*(dirac(omega+E_1+E_2)-dirac(omega-E_1-E_2));
    %intresa = quad2d(resa,-Inf,Inf,-Inf,Inf);
    intresa = vpaintegral(vpaintegral(resa(epsilon_1,epsilon_2),epsilon_1,-Inf,Inf),epsilon_2,-Inf,Inf);
    resb = @(epsilon_1,epsilon_2) (gamma_1-gamma_2).*(1+Delta_inf./(E_1.*E_2)).*(dirac(omega+E_1-E_2)-dirac(omega-E_1+E_2));
    %intresb = quad2d(resb,-Inf,Inf,-Inf,Inf);
    intresb = vpaintegral(vpaintegral(resb(epsilon_1,epsilon_2),epsilon_1,-Inf,Inf),epsilon_2,-Inf,Inf);
    Realsigmaa = 1./(8.*omega).*intresa;
    Realsigmab = 1./(8.*omega).*intresb;
    xaxis(K) = omega./Delta_inf;
    Realsigma(K) = Realsigmaa+Realsigmab;
    plot(xaxis, Realsigma);
    drawnow();
end

[Rebecca Müller]

ok, I found a completely different way to solve my problem. Anyways, thank you really so so much for your help - although I didn't come to the solution that way I feel I learned something. Thanks!