Operator '+' is not supported for operands of type 'function_handle'.

Trying to run the following program

% Define the domain and grid size 
L = 30; % length of the domain 
N = 20; % number of grid points in each direction 
x = linspace(-L/2,L/2,N); 
y = linspace(-L/2,L/2,N); 
z = linspace(-L/2,L/2,N); 
[X,Y,Z] = meshgrid(x,y,z); 
dx = x(2)-x(1); 

% Define the potential function 
V = -1./sqrt(X.^2+Y.^2+Z.^2); % Coulomb potential 

% Define the Laplacian operator using central differences 
Lap = @(f)(-0.5)*(circshift(f,[0,0,-1])+circshift(f,[0,0,1])+circshift(f,[0,-1,0])+circshift(f,[0,1,0])+circshift(f,[-1,0,0])+circshift(f,[1,0,0])-6*f)/(dx^2); 
% Set up the initial wave function (1s state) 
psi = exp(-sqrt(X.^2+Y.^2+Z.^2)); % Gaussian wave packet 

% Normalize the wave function 
psi = psi./sqrt(sum(abs(psi(:)).^2*dx^3)); 

% Define the Hamiltonian operator 
%H=Lap+diag(V(:));
H=Lap+diag(V(:));
%H=@(f)Lap(f)+diag(V(:));
% Set up the time evolution parameters 
dt = 0.01; % time step 
tmax = 10; % maximum time 
t = 0:dt:tmax; 

% Solve the time-dependent Schrodinger equation using the Crank-Nicolson method 
for n = 1:length(t) 
    psi = (eye(N^3)-0.5i*dt*H)\(eye(N^3)+0.5i*dt*H)*psi; 
    psi = psi./sqrt(sum(abs(psi(:)).^2*dx^3)); % renormalize 
end 

% Calculate the energy levels by finding the eigenvalues of the Hamiltonian 
[E,D] = eigs(H,10,'sr'); 
E = diag(E); % extract the eigenvalues 
E = E(E<0); % take only the negative energy levels (bound states) 

% Print the energy levels 
disp('Energy levels:'); 
disp(E); 

I have the follwing

Operator '+' is not supported for operands of type 'function_handle'.

I do not understand how to handle these diabolic objects function_handles.

I have tried all different suggestions I found in the net but simply nothing worked.

Is there something wrong with this code? I copied this from a site.

Thank you

Demetris Ghikas

Answer 1

[Stephen23]

"Is there something wrong with this code?"

You are calling PLUS on the function handle LAP.

Function handles are not numbers, they cannot be added like this:

H = Lap+diag(V(:));

Curiously the code includes a commented-out line which would work:

H = @(f)Lap(f)+diag(V(:));

But after fixing that your code will throw another error due to calling the function H without any input argument:

psi = (eye(N^3)-0.5i*dt*H)\(eye(N^3)+0.5i*dt*H)*psi;

What value of f do you expect the function H to use when you call it like that?

Answer 2

[Kanishk]

Hi @Demetris,

It seems you're encountering an issue with function handles in MATLAB. A function handle is a versatile data type that allows you to reference a function indirectly. This capability is particularly useful for passing functions as arguments, storing them in variables, or defining anonymous functions.

If your goal is to combine the outputs of two functions, note that directly adding function handles is not meaningful. Instead, you can create an anonymous function to achieve this.

H = @(x) Lap(x) + diag(V(:)); 

This approach lets you define a new function handle ‘H’ that computes the sum of the outputs from ‘Lap(x)’ and ‘diag(V(:))’.

Looking at your workflow, you might require the Hamiltonian ‘H’ as a sparse matrix instead of a function handle before using the Crank-Nicolson method in the loop and calculating eigen values. This also requires calculating Laplace ‘Lap’ as a matrix instead of function handle.

To learn more about function handles, please go through this official MATLAB documentation.

https://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html

Hope this helps!

Thanks