# Heaviside function is not integrated

3 views (last 30 days)
Ivo Hoffmanns on 28 Nov 2020
Commented: Ivo Hoffmanns on 28 Nov 2020
The following piece of code does not evaluate the integral.
clear all
close all
%% Castigliano
E=210e9; %e-modulus aluminium=69 steel-210
rho=8050; %density of aluminium kg/m^3=2702, steel=8050
g=9.81;
length=[10,10,10]; %length of the beams
width=[0.5,0.45,0.4]; %width of the 1st, 2nd and 3rd beam
w_thickness=[0.05,0.05,0.05]; %thickness of the vertical sides
h_thickness=[0.05,0.05,0.05]; %thickness of the horizontal sides
height=[0.5,0.45,0.4];
syms A B C x ; %A B and C are loads located where two sections meet and the end
% A=0;
% B=0;
% C=30*9.81;
I1=(width(1)*height(1)^3)-((width(1)-w_thickness(1))*(height(1)-h_thickness(1))^3)/(12); %Moment of inertia
I2=(width(2)*height(2)^3)-((width(2)-w_thickness(2))*(height(2)-h_thickness(2))^3)/(12);
I3=(width(3)*height(3)^3)-((width(3)-w_thickness(3))*(height(3)-h_thickness(3))^3)/(12);
I=heaviside(x)*I1-heaviside(x-length(1))*(I1-I2)-heaviside(x-(length(1)+length(2)))*(I2-I3);
M=A*(length(1)-x)+B*(length(1)+length(2)-x)+C*(length(1)+length(2)+length(3)-x);
Mda=length(1)-x;
Mdb=length(1)+length(2)-x;
Mdc=length(1)+length(2)+length(3)-x;
function3=M*Mdc/(E*I);
delta3=int(function3,x,0,sum(length))
%eval(delta3)
Instead of solving the integral, it gives the following result
delta3 =
int(-((x - 30)*(A*(x - 10) + B*(x - 20) + C*(x - 30)))/((37331305981018095703125*heaviside(x - 10))/8796093022208 + (214619172183736318359375*heaviside(x - 20))/70368744177664 - (436546248401485810546875*heaviside(x))/35184372088832), x, 0, 30)

Alan Stevens on 28 Nov 2020
You could do it this way, without the symbolic stuff:
%% Castigliano
E=210e9; %e-modulus aluminium=69 steel-210
rho=8050; %density of aluminium kg/m^3=2702, steel=8050
g=9.81;
length=[10,10,10]; %length of the beams
width=[0.5,0.45,0.4]; %width of the 1st, 2nd and 3rd beam
w_thickness=[0.05,0.05,0.05]; %thickness of the vertical sides
h_thickness=[0.05,0.05,0.05]; %thickness of the horizontal sides
height=[0.5,0.45,0.4];
%A B and C are loads located where two sections meet and the end
A=0;
B=0;
C=30*9.81;
I1=(width(1)*height(1)^3)-((width(1)-w_thickness(1))*(height(1)-h_thickness(1))^3)/(12); %Moment of inertia
I2=(width(2)*height(2)^3)-((width(2)-w_thickness(2))*(height(2)-h_thickness(2))^3)/(12);
I3=(width(3)*height(3)^3)-((width(3)-w_thickness(3))*(height(3)-h_thickness(3))^3)/(12);
I=@(x) heaviside(x)*I1-heaviside(x-length(1))*(I1-I2)-heaviside(x-(length(1)+length(2)))*(I2-I3);
M=@(x) A*(length(1)-x)+B*(length(1)+length(2)-x)+C*(length(1)+length(2)+length(3)-x);
Mda=@(x) length(1)-x;
Mdb=@(x) length(1)+length(2)-x;
Mdc=@(x) length(1)+length(2)+length(3)-x;
function3=@(x) M(x).*Mdc(x)./(E*I(x));
delta3=integral(function3,0,sum(length));
disp(delta3)

#### 1 Comment

Ivo Hoffmanns on 28 Nov 2020
That works, thanks!

R2020a

### Community Treasure Hunt

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

Start Hunting!

Translated by