Matlab: function handle integration with several variables

My goal here is to build an array (cell array since I'm working with function handles) via a for loop and take the integral of each element, plug in a value and get an array. But I get the following error:

Input function must return 'double' or 'single' values. Found 'function_handle'.

The error occurs on the line when I'm trying to plug in the value 1 (or any scalar value) for x_2. Any tips on how to "handle" this error? Note a(1,1) and c(1,1) are both scalar values (eg. 0 and 1).

Here is the code:

    FUN_1 = @(y_1,y_2,x_1,x_2)sum(heaviside(y_1-a_k(1:m,1)).*dirac(1,y_2-a_k(1:m,2))).*(-1/2.*log((x_1-y_1).^2+(x_2-y_2).^2))+(x_1-y_1).^2./((x_1-y_1).^2)+sum(dirac(y_1-a_k(1:m,1)).*dirac(y_2-a_k(1:m,2))).*(-1/2.*log((x_1-y_1).^2+(x_2-y_2).^2))+(x_1-y_1).*(x_2-y_2)./((x_1-y_1).^2+(x_2-y_2).^2);
Q_1 = @(x_1,x_2)integral2(@(y_1,y_2)FUN_1(y_1,y_2,x_1,x_2),a(1,1),c(1,1),a(1,2),c(1,2));

 FUN_2 = @(y_1,y_2,x_1,x_2)sum(heaviside(y_1-a_k(1:m,1)).*dirac(1,y_2-a_k(1:m,2))).*(-1/2.*log((x_1-y_1).^2+(x_2-y_2).^2))+(x_1-y_1).*(x_2-y_2)./((x_1-y_1).^2)+sum(dirac(y_1-a_k(1:m,1)).*dirac(y_2-a_k(1:m,2))).*(-1/2.*log((x_1-y_1).^2+(x_2-y_2).^2))+(x_2-y_2).^2./((x_1-y_1).^2+(x_2-y_2).^2);
Q_2 = @(x_1,x_2)integral2(@(y_1,y_2)FUN_1(y_1,y_2,x_1,x_2),a(1,1),c(1,1),a(1,2),c(1,2));

    k = zeros(1,2*M);
    n=0;
for n = 0:2*M-1
    S = @(x_1,x_2)Q_1(x_1,x_2)*2*n*(x_1+1i*x_2)^(n-1) + Q_2(x_1,x_2)*2*n*1i*(x_1+1i*x_2)^(n-1);
    R = @(x_2)integral(@(x_1)S,a(1,1),c(1,1));
    k(1,n+1) = R(1);
end
disp(k);
  end

Answer 1

[Guillaume]

A guess as I've not really tried to understand what your functions are doing:

R = @(x_2) integal(@(x_1) S(x_1, x_2), a, c)

The error you get makes sense. (@(x_1) S is a function that returns the function handle S regardless of the input. My modification returns the result of S(x_1, x_2) for input x_1 received from integral.

Note that if a and c are scalar, there's absolutely no point in writing them as a(1, 1) and c(1, 1) other than puzzling the reader and making them wonder if you meant to pass instead some other variable that was a 2D array.

[Pro B]

Right, I made the modification, however after making this modification I get an error in matrix dimension for this line: k =( 1,n+1) = R(1). Would you have any idea as to why? a and c are both 1 by 2 arrays that have to be inputed by the user at the beginning (this is only a part of my program).

[Guillaume]

Well, what is size of

R(1)

when the error occurs. It clearly not scalar from the error message.

[Pro B]

The size is a 1 by 1 function handle.

[Pro B]

Well I tried to simplify my code by making n=1, however I am still getting the same error in matrix dimension when I plug in for R(1). Any ideas?

prompt = 'Enter a st. [x_a y_a]'; % Enter boundaries of domain (lower left hand corner), sze(1,2)
a = input(prompt);
prompt = 'Enter c st. [x_c y_c]'; % Upper right hand corner of the domain, sze(1,2)
c = input(prompt);
prompt = 'Enter m'; % Number of dislocations, scalar
m = input(prompt);
prompt = 'Enter coordinates of dislocations st. [x_1 y_1; ... ; x_m, y_m]'; % sz(m,2)
a_k = input(prompt);

% Computation of the convolution of Phi and F 

FUN_1 = @(y_1,y_2,x_1,x_2)sum(heaviside(y_1-a_k(1:m,1)).*dirac(1,y_2-a_k(1:m,2))).*(-1/2*log((x_1-y_1).^2+(x_2-y_2).^2))+sum(dirac(y_1-a_k(1:m,1)).*dirac(y_2-a_k(1:m,2)));
Q_1 = @(x_1,x_2)integral2(@(y_1,y_2)FUN_1(y_1,y_2,x_1,x_2),a(1,1),c(1,1),a(1,2),c(1,2));
FUN_2 = @(y_1,y_2,x_1,x_2)sum(heaviside(y_1-a_k(1:m,1)).*dirac(1,y_2-a_k(1:m,2)))+sum(dirac(y_1-a_k(1:m,1)).*dirac(y_2-a_k(1:m,2)))*(-1/2*log((x_1-y_1).^2+(x_2-y_2).^2));
Q_2 = @(x_1,x_2)integral2(@(y_1,y_2)FUN_1(y_1,y_2,x_1,x_2),a(1,1),c(1,1),a(1,2),c(1,2));

S = @(x_1,x_2)Q_1(x_1,x_2)*2.*(x_1+1i*x_2)+Q_2(x_1,x_2)*2*1i.*(x_1+1i*x_2);
R = @(x_2)integral(@(x_1)S(x_1,x_2),a(1,1),c(1,1));
b = R(1);
disp(b)

[Torsten]

Aside from the formal MATLAB error from above: it makes little sense to apply numerical integration methods to functions that are full of discontinuities (dirac, heaviside).

Did you try to plot abs((R)) ? How does it look ? Not smooth, I guess ...

Best wishes

Torsten.

[Walter Roberson]

You posted a Question about this; I replied there. You missed that integral() passes in a vector of values.