Problems with integration using functions

Dear all, I have a problem when doing integration twice with a defined function. My code is as follows: First I define a function psitimesX(x) with variable x, this function has to be integrated over variable p by a function handle under function psitimesX. Next I call back the function and integrate it over x from -infinity to +infinity. I am not sure what has does wrongly as matlab just does not gives me input(need to wait for a long time), I think there is something wrong. Or is there an alternative way to do this task of integrating twice?

clear all;
%%%integrating variable over x
expecX_xbasis = integral(@psitimesX,-inf, inf,'ArrayValued',true)

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

function psitimesX = psitimesX(x)
     %%%Define parameters
     J = -0.4 ;
     theta = 0.4;
     omega_y =0.04;
     omega_z = -0.052;
     m=1/(-2*J);
     gammavar= 5;
     gamma_profile=0.05;
     correctionfactor = (pi/gamma_profile)^0.25;
     %%%function of x, with function handle with variable p
     wavefunction_xbasis_component1 = @(p)(1/sqrt(2)).*(1./sqrt(2.*pi)).*correctionfactor.*((exp(-1i.*(0:0.1:150).*(2.*J*cos(p).*cos(theta))-gammavar.*p.^2 + 1i.*p.*x))).*(cos((0:0.1:150).*sqrt(omega_y^2 +omega_z^2 + J^2 -(J^2).*cos(2*p) -2.*J.*sin(p)*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta)))  ) + 1i.*(omega_z - 2.*J*sin(p).*sin(theta)).*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2)*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))/sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta)*sin(p) + 2.*omega_z*sin(theta))) +omega_y.*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))./sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))    ;
     wavefunction_xbasis_component2 = @(p)(1/sqrt(2)).*(1./sqrt(2.*pi)).*correctionfactor.*((exp(-1i.*(0:0.1:150).*(2.*J.*cos(p).*cos(theta))-gammavar.*p.^2 + 1i.*p.*x))).*(cos((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2*theta).*sin(p) + 2.*omega_z.*sin(theta)))) -1i.*(omega_z - 2.*J.*sin(p).*sin(theta)).*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p)*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))./sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p)*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))) -omega_y.*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta)*sin(p) + 2.*omega_z.*sin(theta))))./sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2*omega_z.*sin(theta))))    ;
     %%%psitimesX is a function of x, after integrating over variable p
     psitimesX = x.*(abs(integral(wavefunction_xbasis_component1,-inf, inf,'ArrayValued',true)).^2 + abs(integral(wavefunction_xbasis_component2,-inf, inf,'ArrayValued',true)).^2 ) ;

end

[Star Strider]

With this slight change:

expecX_xbasis = integral(@(x)psitimesX(x),-inf, inf,'ArrayValued',true)

I was able to get it to work, however it is taking forever and throwing this warning:

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.5e-06. The

integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

three times in the first 15 miinutes of running time.

I am stopping it.

You can run it and then see what the problems are wiith it.

clear all;
%%%integrating variable over x
expecX_xbasis = integral(@(x)psitimesX(x),-inf, inf,'ArrayValued',true)

function psitimesX = psitimesX(x)
     %%%Define parameters
     J = -0.4 ;
     theta = 0.4;
     omega_y =0.04;
     omega_z = -0.052;
     m=1/(-2*J);
     gammavar= 5;
     gamma_profile=0.05;
     correctionfactor = (pi/gamma_profile)^0.25;
     %%%function of x, with function handle with variable p
     wavefunction_xbasis_component1 = @(p)(1/sqrt(2)).*(1./sqrt(2.*pi)).*correctionfactor.*((exp(-1i.*(0:0.1:150).*(2.*J*cos(p).*cos(theta))-gammavar.*p.^2 + 1i.*p.*x))).*(cos((0:0.1:150).*sqrt(omega_y^2 +omega_z^2 + J^2 -(J^2).*cos(2*p) -2.*J.*sin(p)*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta)))  ) + 1i.*(omega_z - 2.*J*sin(p).*sin(theta)).*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2)*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))/sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta)*sin(p) + 2.*omega_z*sin(theta))) +omega_y.*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))./sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))    ;
     wavefunction_xbasis_component2 = @(p)(1/sqrt(2)).*(1./sqrt(2.*pi)).*correctionfactor.*((exp(-1i.*(0:0.1:150).*(2.*J.*cos(p).*cos(theta))-gammavar.*p.^2 + 1i.*p.*x))).*(cos((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2*theta).*sin(p) + 2.*omega_z.*sin(theta)))) -1i.*(omega_z - 2.*J.*sin(p).*sin(theta)).*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p)*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))))./sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p)*(J.*cos(2.*theta).*sin(p) + 2.*omega_z.*sin(theta))) -omega_y.*sin((0:0.1:150).*sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta)*sin(p) + 2.*omega_z.*sin(theta))))./sqrt(omega_y.^2 +omega_z.^2 + J.^2 -(J.^2).*cos(2.*p) -2.*J.*sin(p).*(J.*cos(2.*theta).*sin(p) + 2*omega_z.*sin(theta))))    ;
     %%%psitimesX is a function of x, after integrating over variable p
     psitimesX = x.*(abs(integral(wavefunction_xbasis_component1,-inf, inf,'ArrayValued',true)).^2 + abs(integral(wavefunction_xbasis_component2,-inf, inf,'ArrayValued',true)).^2 ) ;
end

.

[Tsz Tsun]

It runs forever, I believe there should be some problems with the code, as I suppose numerical integration can run effectively? I am thinking if there is another way to write it.

[Star Strider]

I am runniing it in the background, on MATLAB Online. It appears to me to be not close to converging.

Answer 1

[Abhishek Kumar Singh]

Hi @Tsz Tsun,

The issue you're encountering seems to arise from a few potential problems in your MATLAB code related to integration and function handling. Here's a breakdown of the problems and how you can address them:

1. Function Definition and Usage: Ensure that the psitimesX function is properly defined and accessible. If you encounter errors related to function calls, such as the one you experienced, consider using an anonymous function to correctly pass variables, as in @(x) psitimesX(x).
2. Infinite Limits in Integration: Using -inf and inf as limits can lead to convergence issues, especially if the function being integrated does not decay to zero at infinity. Consider using finite limits if possible or ensure the function is well-behaved over infinite limits.
3. Array Valued Integration: The 'ArrayValued', true option is used when the integrand returns an array for each input. Ensure that your function handle returns a scalar for each input p and x if this option is not intended.
4. Complex Integrands: If your integrand is complex, ensure that the integration process can handle complex numbers. MATLAB's integral function can handle complex values, but make sure the integrand is properly defined.
5. Vectorization: Ensure that the function handles are vectorized properly. MATLAB's integral function is designed to work with vectorized functions, which means your function should be able to accept and return vectors.

Here's how you can structure your code to address these issues:

% Define the function handle for psitimesX
psitimesX_handle = @(x) psitimesX(x);

% Perform the integration over x
expecX_xbasis = integral(psitimesX_handle, -10, 10, 'ArrayValued', true);  % Use finite limits for testing

% Function definition
function result = psitimesX(x)
    % Define parameters
    J = -0.4;
    theta = 0.4;
    omega_y = 0.04;
    omega_z = -0.052;
    m = 1 / (-2 * J);
    gammavar = 5;
    gamma_profile = 0.05;
    correctionfactor = (pi / gamma_profile)^0.25;

    % Define function handles for wavefunction components
    wavefunction_xbasis_component1 = @(p) (1 / sqrt(2)) * (1 / sqrt(2 * pi)) * correctionfactor .* ...
        exp(-1i * (0:0.1:150) .* (2 * J * cos(p) * cos(theta)) - gammavar * p.^2 + 1i * p * x) .* ...
        (cos((0:0.1:150) .* sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta)))) + ...
        1i * (omega_z - 2 * J * sin(p) * sin(theta)) .* sin((0:0.1:150) .* sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta)))) ./ ...
        sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta))) + ...
        omega_y * sin((0:0.1:150) .* sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta)))) ./ ...
        sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta))));

    wavefunction_xbasis_component2 = @(p) (1 / sqrt(2)) * (1 / sqrt(2 * pi)) * correctionfactor .* ...
        exp(-1i * (0:0.1:150) .* (2 * J * cos(p) * cos(theta)) - gammavar * p.^2 + 1i * p * x) .* ...
        (cos((0:0.1:150) .* sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta)))) - ...
        1i * (omega_z - 2 * J * sin(p) * sin(theta)) .* sin((0:0.1:150) .* sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta)))) ./ ...
        sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta))) - ...
        omega_y * sin((0:0.1:150) .* sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta)))) ./ ...
        sqrt(omega_y^2 + omega_z^2 + J^2 - (J^2) * cos(2 * p) - 2 * J * sin(p) * (J * cos(2 * theta) * sin(p) + 2 * omega_z * sin(theta))));

    % Integrate over p
    integral1 = integral(wavefunction_xbasis_component1, -inf, inf, 'ArrayValued', true);
    integral2 = integral(wavefunction_xbasis_component2, -inf, inf, 'ArrayValued', true);

    % Calculate result
    result = x .* (abs(integral1).^2 + abs(integral2).^2);
end

Checking the final output:

plot(expecX_xbasis)

I hope this helps!

[Tsz Tsun]

@Abhishek Kumar Singh Thanks a lot! It helps!

[Walter Roberson]

MATLAB requires functions to be defined in their own files or as nested functions.

This is incorrect. It is fine to define functions at the end of function files: they will be treated as private functions. They do not need to be defined in their own files or as nested functions.

Recently (R2024a I think) it also became file to define functions in the middle of script files, provided that you are not within any kind of control structure.

[Abhishek Kumar Singh]

Thank @Walter Roberson for the correction- I have edited that part now. I checked documentation and indeed from R2016b we can define functions at end of scripts as local functions and can be added anywhere in the file except within conditional contexts. Before R2024a there was a constraint that local functions in scripts must be defined at the end of the file, after the last line of script code.