numerical integration 2d function of three variables and plotting

Question

Rebecca Müller 2018-6-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/407092-numerical-integration-2d-function-of-three-variables-and-plotting

评论： Rebecca Müller 2018-6-24

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);

2 个评论
显示无隐藏无

Jan 2018-6-23

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

Walter Roberson 2018-6-23

quad2d cannot be used for infinite ranges.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2018-6-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/407092-numerical-integration-2d-function-of-three-variables-and-plotting#answer_325996

在 MATLAB Online 中打开

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);

8 个评论
显示 6更早的评论隐藏 6更早的评论

Walter Roberson 2018-6-24

在 MATLAB Online 中打开

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 2018-6-24

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!

请先登录，再进行评论。