Function handle integration with several variable, matrix dimension error

Hi! I am using function handles to integrate with several variables (note: I made the same program using symbolic integration, but it takes a really long time). However, I keep getting an error in matrix dimensions in the last line of my when I plug in R(1) and I really can't figure out what is wrong. Any tips would be greatly appreciated!

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

FUN_2 = @(y_1,y_2,x_1,x_2)sum(heaviside(y_1-1).*dirac(1,y_2-1))+sum(dirac(y_1-1).*dirac(y_2-1))*(-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),1,2,1,2);

S = @(x_1,x_2)Q_1(x_1,x_2)+Q_2(x_1,x_2);
R = @(x_2)integral(@(x_1)S(x_1,x_2),1,2);
b = R(1);
disp(b)

Answer 1

[Walter Roberson]

integral() calls the function passing in a vector of values. So your S will be called with a vector x_1, so your Q_1 will be called with a vector x_1. Your Q_1 calls integral2, which is going to call FUN_1 with array y_1 and y_2 (the same size) and the vector x_1 (of different size). Your FUN_1 expects x_1 and y_1 to be compatible size for subtraction, but the two are not compatible.

You could use the 'ArrayValued' option of integral so that it only passes in a scalar for x_1, at a cost to efficiency. You should look more carefully at which operations can be vectorized. Sometimes it helps to use arrayfun() along the way.

[Pro B]

Thank you. I now understand where the error comes from and I was able to solve it by using the 'ArrayValued' on x_1 so that it would of compatible size with y_1. My code works now, but the results of the integrals I get are wrong. For instance Q_1 and Q_2 always give me 0 if I plug in any random value eg. x_1 = 0.4 and x_2 = 0.4. (I know for a fact that this answer is wrong since on the program where I do symbolic integration I obtain Q_1 = atan(8/3) - pi/4 and Q_2 = -log(18/25)/2. Is this error do the fact that I am not using function handles properly?

[Pro B]

This is the code with symbolic expressions:

syms x_1 x_2 y_1 y_2
F_1(y_1, y_2) = sum(heaviside(y_1-1).*dirac(1,y_2-1));
F_2(y_1,y_2) = sum(dirac(y_1-1).*dirac(y_2-1));
Phi(x_1,x_2) = -1/2*log((x_1-y_1).^2 + (x_2-y_2).^2)); % Phi represents the fundamental tensor

% Computation of the convolution of Phi and F

w_1 = int(F_1(y_1,y_2).*Phi(x_1,x_2), y_2);
u_1 = subs(w_1,y_2, 2) - subs(w_1, y_2, 0);
w_1 = int(u_1, y_1);
w_1 = simplify(w_1);
v_1 = subs(w_1, y_1, 2) - subs(w_1, y_1, 0);

w_2 = int(F_2(y_1,y_2)*Phi(x_1,x_2), y_2);
u_2 = subs(w_2,y_2, 2) - subs(w_2, y_2, 0);
w_2 = int(u_2, y_1);
w_2 = simplify(w_2);
v_2 = subs(w_2, y_1, 2) - subs(w_2, y_1, 0);

V = [v_1; v_2];
disp(V)
vars = [x_1,x_2];
newVars = [0.4,0.4];
Y = subs(V, vars, newVars);
disp(Y) %Y represents the same thing as Q_1 and Q_2

result: Y= [atan(8/3) - pi/4; -log(18/25)/2]

This is what I get with my new code:

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

FUN_2 = @(y_1,y_2,x_1,x_2)sum(dirac(y_1-1).*dirac(y_2-1)).*(-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),0,2,0,2);

Q_1(0.4,0.4)
Q_2(0.4,0.4)

Q_1 = 0 and Q_2 = 0

[Walter Roberson]

Your dirac(y_2-1) and dirac(1,y_2-1) will only ever be non-zero if y_2 == 1 exactly, which is not at all certain to occur when you are doing 2D integration.

When you use integral() calls, but not when you use integral2() calls, you can use the Waypoints option to force evaluation at particular points.

But if you did manage to force evaluation at y_2 == 1 exactly then the result of the dirac() would be inf, and the result of the dirac(1) would be -inf . Those would then be multiplied by the heaviside(y_1-1), which would be 0 for y_1 < 1 -- and inf or -inf times 0 is NaN. Having a NaN in a sum "pollutes" the sum, forcing it to be NaN, so you would end up with a NaN result.

You should, in short, never do numeric integration with a dirac() term. You also need to worry about doing numeric integration with a heaviside term because the formal definition of Heaviside leaves the value at 0 undefined -- so you need to know which value heaviside(0) is going to return. For the last couple of releases, the value of heaviside(0) can be configured by using sympref(); the default is 1/2 -- but it is common for people to use heaviside when they mean the unit step function, where u(0) has been defined as 1.

What you should do instead of numeric integration with a dirac term (and ideally, with a heaviside term too) is break the integration up into as many pieces as needed, and add the results together, taking care to handle the exact boundaries in a manner that is meaningful for your situation.

[Pro B]

Yes, after posting my question I quickly realized what was wrong with numeric integration for the dirac and heaviside functions. I was able to eventually finish the program, by performing the integrations involving the Dirac delta and Heaviside function by hand and putting the result of the integration directly into the program. Thank you very much for all your help!