ODE45 returns NaN values.

Question

Bhanu Pratap Akherya 2021-8-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1439224-ode45-returns-nan-values

评论： Star Strider 2021-8-24

l&d.xlsx

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2021-8-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1439224-ode45-returns-nan-values#answer_773059

在 MATLAB Online 中打开

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.

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2021-8-23

在 MATLAB Online 中打开

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 2021-8-24

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 2021-8-24

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’.

.

请先登录，再进行评论。