Adding:
% Qt = [t>=min(t_array) t<=max(t_array)]
to the ‘beam_function’ code demonstrates the problem. The ‘t’ value is always greater than the highest value of ‘t_array’ so interp1 returns NaN since it is not instructed on how to extrapolate. Adding that capability, and changing the solver to ode15s (since this is apparently a ‘stiff’ system) returns these results —
% global K M C u;
Ne=6;
l=1; %length
t=0.02; %thickness
b=0.02; %width
modulus=2e11; %(E)
area=b*t;
imoment=(b*((t)^3))/12;
Le=l/Ne; %length of element
Rho=7850; %density
%Element stiffness matrix
K1=(modulus*imoment/(Le^3))*[12,6*Le,-12,6*Le; ...
6*Le,4*Le*Le,-6*Le,2*Le*Le; ...
-12,-6*Le,12,-6*Le; ...
6*Le,2*Le*Le,-6*Le,4*Le*Le];
Kglobal=zeros(2*(Ne+1),2*(Ne+1));
M1=[156 22*Le 54 -13*Le;...
22*Le 4*Le*Le 13*Le -3*Le*Le;...
54 13*Le 156 -22*Le;...
-13*Le -3*Le*Le -22*Le 4*Le*Le]*(Rho*Le*b*t)/420;
Mglobal=zeros(2*(Ne+1),2*(Ne+1));
for ii=1:Ne
Kglobal(2*ii-1:2*(ii+1),2*ii-1:2*(ii+1))=Kglobal(2*ii-1:2*(ii+1),2*ii-1:2*(ii+1))+K1;
Mglobal(2*ii-1:2*(ii+1),2*ii-1:2*(ii+1))=Mglobal(2*ii-1:2*(ii+1),2*ii-1:2*(ii+1))+M1;
end
K=Kglobal;
K(1:2,:)=[];
K(:,1:2)=[];
M=Mglobal;
M(1:2,:)=[];
M(:,1:2)=[];
C=0.05*Kglobal;
C(1:2,:)=[];
C(:,1:2)=[];
K
M
C
u=(2*Ne)+1;
dt=0.001;
T=300;
%Displacement initials
y0=zeros(2*(2*(Ne+1))-4,1);
y0(end-1,1)=0.5;
%ODE function
% a=xlsread('l&d.xlsx');
a = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/719254/l&d.xlsx')
t_array = a(1,:); % This is t array from xls file
f_array = a(2,:); % This is F array from xls file
[tsol ysol]=ode15s(@(t, y) beam_function(t, y, t_array, f_array, K, M, C, u),[1:dt:T],y0);
plot(tsol,ysol(:,Ne))
function [dy]=beam_function(t,y, t_array, f_array, K, M, C, u)
% global K M C u;
F = interp1(t_array,f_array,t, 'linear','extrap');
dy=[y(u:end);
M\(F-K*y(1:u-1)-C*y(u:end))];
end
If you want different results, it will be necessary to scale ‘t’ to be within the limits of ‘t_array’ so that the interpolation works without the need to extrapolate.
I also eliminated the global variables and passed them as extra parameters to ‘beam_funciton’. See Passing Extra Parameters for details.
.
