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

3 次查看（过去 30 天）
Diego Abad
Diego Abad 2024-5-29
编辑： Dyuman Joshi 2024-5-29
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.

请先登录，再回答此问题。

采纳的回答

Dyuman Joshi
Dyuman Joshi 2024-5-29
编辑：Dyuman Joshi 2024-5-29
Ran in:
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

更多回答（1 个）

the cyclist
the cyclist 2024-5-29
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.
  2 个评论
the cyclist
the cyclist 2024-5-29
Ran in:
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
the cyclist 2024-5-29
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Install Products 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by