ODE45 returns NaN values.

The ODE45 function is returning a NaN value for the dy. I am a beginner at MATLAB coding, I do not know where the issue is. Can anyone help? Also I have attached the .xlsx file.

Here is my code:

Main code:

clear all

clc

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');

t_array = a(1,:); % This is t array from xls file

f_array = a(2,:); % This is F array from xls file

[tsol ysol]=ode45(@(t, y) beam_function(t, y, t_array, f_array),[1:dt:T],y0);

plot(tsol,ysol(:,Ne))

Function code:

function [dy]=beam_function(t,y, t_array, f_array)

global K M C u;

F = interp1(t_array,f_array,t);

dy=[y(u:end);

M\(F-K*y(1:u-1)-C*y(u:end))]

Answer 1

[Star Strider]

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.

.

[Bhanu Pratap Akherya]

thanks for your help, can you also show the results you got here?

[Star Strider]

As always, my pleasure!

I did originally, however I corrected some errors and forgot to re-run the code, so I’ll so it here —

% 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

K = 12×12

1.0e+07 * 1.3824 0 -0.6912 0.0576 0 0 0 0 0 0 0 0 0 0.0128 -0.0576 0.0032 0 0 0 0 0 0 0 0 -0.6912 -0.0576 1.3824 0 -0.6912 0.0576 0 0 0 0 0 0 0.0576 0.0032 0 0.0128 -0.0576 0.0032 0 0 0 0 0 0 0 0 -0.6912 -0.0576 1.3824 0 -0.6912 0.0576 0 0 0 0 0 0 0.0576 0.0032 0 0.0128 -0.0576 0.0032 0 0 0 0 0 0 0 0 -0.6912 -0.0576 1.3824 0 -0.6912 0.0576 0 0 0 0 0 0 0.0576 0.0032 0 0.0128 -0.0576 0.0032 0 0 0 0 0 0 0 0 -0.6912 -0.0576 1.3824 0 -0.6912 0.0576 0 0 0 0 0 0 0.0576 0.0032 0 0.0128 -0.0576 0.0032

M

M = 12×12

0.3888 0 0.0673 -0.0027 0 0 0 0 0 0 0 0 0 0.0003 0.0027 -0.0001 0 0 0 0 0 0 0 0 0.0673 0.0027 0.3888 0 0.0673 -0.0027 0 0 0 0 0 0 -0.0027 -0.0001 0 0.0003 0.0027 -0.0001 0 0 0 0 0 0 0 0 0.0673 0.0027 0.3888 0 0.0673 -0.0027 0 0 0 0 0 0 -0.0027 -0.0001 0 0.0003 0.0027 -0.0001 0 0 0 0 0 0 0 0 0.0673 0.0027 0.3888 0 0.0673 -0.0027 0 0 0 0 0 0 -0.0027 -0.0001 0 0.0003 0.0027 -0.0001 0 0 0 0 0 0 0 0 0.0673 0.0027 0.3888 0 0.0673 -0.0027 0 0 0 0 0 0 -0.0027 -0.0001 0 0.0003 0.0027 -0.0001

C

C = 12×12

1.0e+05 * 6.9120 0 -3.4560 0.2880 0 0 0 0 0 0 0 0 0 0.0640 -0.2880 0.0160 0 0 0 0 0 0 0 0 -3.4560 -0.2880 6.9120 0 -3.4560 0.2880 0 0 0 0 0 0 0.2880 0.0160 0 0.0640 -0.2880 0.0160 0 0 0 0 0 0 0 0 -3.4560 -0.2880 6.9120 0 -3.4560 0.2880 0 0 0 0 0 0 0.2880 0.0160 0 0.0640 -0.2880 0.0160 0 0 0 0 0 0 0 0 -3.4560 -0.2880 6.9120 0 -3.4560 0.2880 0 0 0 0 0 0 0.2880 0.0160 0 0.0640 -0.2880 0.0160 0 0 0 0 0 0 0 0 -3.4560 -0.2880 6.9120 0 -3.4560 0.2880 0 0 0 0 0 0 0.2880 0.0160 0 0.0640 -0.2880 0.0160

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')

a = 2×300

0.0010 0.0020 0.0030 0.0040 0.0050 0.0060 0.0070 0.0080 0.0090 0.0100 0.0110 0.0120 0.0130 0.0140 0.0150 0.0160 0.0170 0.0180 0.0190 0.0200 0.0210 0.0220 0.0230 0.0240 0.0250 0.0260 0.0270 0.0280 0.0290 0.0300 5.5347 1.3163 1.6579 1.9464 2.2100 2.2922 2.2443 2.1359 2.0355 1.9556 1.8888 1.8303 1.7822 1.7422 1.7064 1.6718 1.6390 1.6105 1.5868 1.5663 1.5475 1.5292 1.5110 1.4927 1.4740 1.4556 1.4375 1.4200 1.4036 1.3888

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

.

[Bhanu Pratap Akherya]

Thanks for your help, I wanted to know one more thing, How does this code work? I mean does F value take 1 value of f_array and t_array interpolates it and feeds it to the dy equation at one ODE time step and then takes the next values of f_array and t_array interpolates it and feeds it to the dy equation at another ODE time step and so on?

[Star Strider]

As always, my pleasure!

Essentially, yes. The ‘F’ value is interpolated (or extrapolated) from the existing data vectors to the current value of ‘t’ passed to it from ode15s in each call to it. So, it returns one value interpolated (or extrapolated) from ‘t_array’ and ‘f_array’ for each value of ‘t’ presented to it in each call to ‘beam_function’.

.