Problems with integration using functions
17 ビュー (過去 30 日間)
古いコメントを表示
Tsz Tsun
2024 年 9 月 5 日 14:26
編集済み: Abhishek Kumar Singh
2024 年 9 月 5 日 22:58
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)
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
3 件のコメント
Star Strider
2024 年 9 月 5 日 15:57
移動済み: Star Strider
2024 年 9 月 5 日 15:57
I am runniing it in the background, on MATLAB Online. It appears to me to be not close to converging.
採用された回答
Abhishek Kumar Singh
2024 年 9 月 5 日 16:56
編集済み: Abhishek Kumar Singh
2024 年 9 月 5 日 22:58
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:
- 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).
- 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.
- 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.
- 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.
- 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!
3 件のコメント
Walter Roberson
2024 年 9 月 5 日 20:51
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
2024 年 9 月 5 日 22:57
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.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Function Creation についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!