Why I am not getting the same result for an integral of a piecewise function?

I set the following piecewise function to integrate it:

function lambda = weight(x,L1,L2,M)
% L1, L2, M are constants
if x < L1
    lambda = (M ./ (4 .* (L1 + L2))) + ((3 .* M) / (2 .* L2 .^ 2)) .* (L1 - x);
else
    lambda = M ./ (4 .* (L1 + L2)) .* ones(size(x));
end
end

Now, I am trying to integrate lambda between 0 and 15. Everything looks okay, but the results are different when I try to integrate the whole interval than when I do it by sections. In the following code, check should be 0 or near to 0, but I am getting 5.1042e+03.

lambda = @(x) weight(x,5,10,35000);
aux1 = integral(lambda,0,15)
aux2 = integral(lambda,0,5)
aux3 = integral(lambda,5,15)
check1 = aux1 - aux2 + aux3

Does anyone know what is the issue here?

I know it work fine by sections but I would like it to do it as a whole.

Answer 1

[Dyuman Joshi]

1 - It should be aux1 - (aux2 + aux3)

2 - The function is not vectorized properly, which provides incorrect results. I have modified the function to include vectorization, see below -

lambda = @(x) weightvec(x,5,10,35000);

aux1 = integral(lambda,0,15)

aux1 = 1.5313e+04

aux2 = integral(lambda,0,5)

aux2 = 9.4792e+03

aux3 = integral(lambda,5,15)

aux3 = 5.8333e+03

check1 = aux1 - (aux2 + aux3)

check1 = 0.0027

This is a numerical answer with relatively low accuracy i.e. the value is near 0. (You can increase the accuracy, but only upto a certain value/limit - Refer to the documentation of integral for more information regarding this)

Let's try it symbolically utilizing the piecewise function -

L1 = 5; L2 = 10; M = 35000;
syms x
y = piecewise(x<L1, (M ./ (4 .* (L1 + L2))) + ((3 .* M) / (2 .* L2 .^ 2)) .* (L1 - x), ...
    M ./ (4 .* (L1 + L2)));

aux1 = int(y,x,0,15);
aux2 = int(y,x,0,5);
aux3 = int(y,x,5,15);
out = aux1 - (aux2+aux3)

out = 

0

Well, as expected.

%Function with vectorization
function lambda = weightvec(x,L1,L2,M)
lambda = zeros(size(x));
idx =  x < L1;
lambda(idx) = (M ./ (4 .* (L1 + L2))) + ((3 .* M) / (2 .* L2 .^ 2)) .* (L1 - x(idx));
lambda(~idx) = M ./ (4 .* (L1 + L2)) .* ones(1,nnz(~idx));
end

Answer 2

[the cyclist]

You have two mistakes in this code. First, you got the signs wrong on the check. It should be

check1 = aux1 - (aux2 + aux3)

Second, the syntax

if x < L1

is probably not doing what you expect. If x is a vector, this is not going to check each element of x against L1, and calculate the corresponding lambda accordingly. It will only go into the "if" portion of the code when all the elements of x are less than L1.

[the cyclist]

The following is not the most efficient, but I wanted to show that looping over x will give the correct result. (I also changed the relative tolerance, to show that you get very close to zero.)

lambda = @(x) weight(x,5,10,35000);
aux1 = integral(lambda,0,15,"RelTol",1e-12)

aux1 = 1.5312e+04

aux2 = integral(lambda,0,5,"RelTol",1e-12)

aux2 = 9.4792e+03

aux3 = integral(lambda,5,15,"RelTol",1e-12)

aux3 = 5.8333e+03

check1 = aux1 - (aux2 + aux3)

check1 = -1.7499e-09

function lambda = weight(x,L1,L2,M)
    lambda = zeros(size(x));
    for n = 1:numel(x)
        if x(n) < L1
            lambda(n) = (M ./ (4 .* (L1 + L2))) + ((3 .* M) ./ (2 .* L2 .^ 2)) .* (L1 - x(n));
        else
            lambda(n) = M ./ (4 .* (L1 + L2));
        end
    end
end

[the cyclist]

Here is a more efficient version of your function:

function lambda = weight(x,L1,L2,M)
    % L1, L2, M are constants
    lambda = (M ./ (4 .* (L1 + L2))) + (x<L1).*((3 .* M) ./ (2 .* L2 .^ 2)) .* (L1 - x);
end