# ODE45 returns NaN values.

4 views (last 30 days)
Bhanu Pratap Akherya on 23 Aug 2021
Commented: Star Strider on 24 Aug 2021
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
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))]

Star Strider on 23 Aug 2021
% 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
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.
.
Star Strider on 24 Aug 2021
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’.
.