changing values of RHS with each time step in ODE

i have a set of values for the force term "F" (in the equation my"+cy'+ky=F) saved in an excel file which i have recalled using:

F=xlsread('l&d.xlsx','O1:O300');
My main code looks something like:
y0=initial conditions
[tsol ysol]=ode15s('beam_function',[1:dt:T],y0);
plot(tsol,ysol);

whereas, my function code looks something like:
function [dy]=beam_function(t,y)
dy=[y(1:u-1);
    inv(M)*(F-K*y-C*y)]

The problem is that i want to recall each value of F at different time steps but i dont know how to do that, can anyone guide me?

Answer 1

[Wan Ji]

Is F time dependent ? If so, just write a function named Force

function F = Force(t)
t_array = []; % This is t array from xls file
f_array = []; % This is F array from xls file
F = interp1(t_array,f_array,t);
end

And the odefun becomes

function [dy]=beam_function(t,y)
dy=[y(1:u-1);
    inv(M)*(Force(t)-K*y-C*y)];
end

[Bhanu Pratap Akherya]

sorry I am just a beginner to MATLAB but you asked me to do this right?

main code:

dt=0.001;

T=300;

[tsol ysol]=ode15s('beam_function',[1:dt:T],y0);

function code:

function [dy]=beam_function(t,y)

global u K M C

function F = Force(t)

t_array = xlsread('l&d.xlsx','Q1:Q300'); % This is t array from xls file

f_array = xlsread('l&d.xlsx','O1:O300'); % This is F array from xls file

F = interp1(t_array,f_array,t);

end

dy=[y(1:u-1);

inv(M)*(Force(t)-K*y-C*y)]

end

[Wan Ji]

You can take a look at my code on dynamic Bar

function dynamicBar
%% This is the implementation of 1D bar (dynamic)
% using finite element method
% 2 nodes for each element
% Exam 1
clc
L = 10; % total sol domain length
A = 0.1; % cross section
emod = 206e9; % Young's Modulus
Nele = 200; % number of element
node = linspace(0,L,Nele+1)'; % node (coordinate)
Nnode = size(node,1); % node number
ndof = 1; % number of degrees of freedom
% each column represents [node1,node2,element length, element cross area, Young's Modulus]
element = [(1:1:Nele)',(2:1:Nele+1)',diff(node),A*ones(Nele,1),emod*ones(Nele,1)];
density = 7850;
%% obtain the total Stiffness
Ktot = zeros(Nnode*ndof, Nnode*ndof);
for i = 1:1:Nele
    inode = element(i,1);
    jnode = element(i,2);
    elementLength = element(i,3);
    elementArea = element(i,4);
    elementEmod= element(i,5);
    Ke = [1,-1;-1,1]*elementArea*elementEmod/elementLength;
    Ktot([inode,jnode],[inode,jnode]) = Ktot([inode,jnode],[inode,jnode]) + Ke;
end
K = sparse(Ktot);
mass = element(:,3).*element(:,4).*density;
massVec = [movsum(mass(:,1)',2)/2, mass(end,1)/2];
M = sparse(1:Nnode*ndof, 1:Nnode*ndof, massVec, Nnode*ndof,Nnode*ndof);
alpha = 0.2;
beta = 0.03;
% alpha = 0;
% beta = 0;
C = alpha*M + beta*K;
F = sparse(Nnode*ndof,1,1,Nnode*ndof,1); % a unit force on the end
dt = 0.00001;
T = 0.01;
N = Nnode*ndof;
y0 = sparse(2*N,1);
ForceFun = @(t) sin(2e4*t)*F; % I dont know your data so I use a function
% You can also use
% Force_from_table = xlsread(...);
% Time_from_table = xlsread(...);
% ForceFun = @(t) interp1(Time_from_table, Force_from_table, t);
[tsol, ysol]=ode23s(@(t,y)barfun(t,y,M,C,K,@(t)ForceFun(t)),0:dt:T, y0);
T = table(tsol, ysol);
save('solution','T');
end

function dydt = barfun(t, y, M, C, K, Force)
N = size(y,1)/2;
y1 = y(1:N,1);
y2 = y(N+1:end,1);
y1(1) = 0; % fixed at x=0
y2(1) = 0;
% x1 = x; x2 = x'
dydt = [y2;
    (Force(t)-K*y1-C*y2)./diag(M)];
% IF M is not diag, Then use following instead
% dydt = [y2;
%    M \ (Force(t)-K*y1-C*y2)];
end

And the second file is for post processing

clc;clear
load('solution.mat');
q = max(max(abs(T.ysol(:,1:end/2))));
N = size(T.ysol,2)/2;
for i = 1:2:size(T,1)
    if(i>1)
    q = (max(abs(T.ysol(i,1:end/2))));
    end
    plot(T.ysol(i,1:end/2))
    axis([0, N, -q, q])
    pause(0.01)
end

All these may help you to some extent.

[Wan Ji]

@Bhanu Pratap Akherya

[Wan Ji]

I have modified your code with following

function main
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)=[];
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; % Here you need make some change!!!!
 
%ODE function
dt2 = 0.1;
t_array = 0:dt2:300; % This is t array from xls file
f_array = sin(2e2*t_array); % This is F array from xls file
[tsol, ysol]=ode15s(@(t, y) beam_function(t, y, t_array, f_array),1:dt:T,y0);
plot(tsol,ysol(:,Ne))
end

function [dy]=beam_function(t,y, t_array, f_array)
global K M C u;
F = interp1(t_array,f_array,t);
% here i have modified for you
dis = y(1:u-1); % displacement
vel = y(u:end); % velocity
dy=[vel;
    M\(F-K*dis-C*vel)];
end

[Bhanu Pratap Akherya]

shouldn't the velocity and displacement terms be interchanged? since the velocity comes first in the initial condition vector?

[Wan Ji]

You know that

[displacement; velocity] = [displacement; velocity] + [velocity; acceleration]*dt

This is the time integration schedule.@Bhanu Pratap Akherya

So here dy is something [velocity; acceleration] here I show.

[Wan Ji]

So, you initial conditions should be modified accordingly.

Or else, if you dont want to change the initial conditions, you need write the beam_function as

function [dy]=beam_function(t,y, t_array, f_array)
global K M C u;
F = interp1(t_array,f_array,t);
% here i have modified for you
vel = y(1:u-1); % velocity
dis = y(u:end); % displacement
dy=[M\(F-K*dis-C*vel);
    vel];
end

All in all, the beam_function you wrote is not true.

[Bhanu Pratap Akherya]

understood... one more thing why did you use that equation for f_array?

[Wan Ji]

Hi, Bhanu Pratap Akherya

Use

F = interp1(t_array,f_array,t);

is only for obtaining the force as time given. t_array and f_array is from your excel, whilst I didnot even see its format.

Answer 2

[Steven Lord]

Use the "ODE with Time-Dependent Terms" example on the documentation page for the ode45 function as a model. The f and g variables in that example correspond to the data from your spreadsheet, with ft and gt being the times associated with that data. Then inside the myode function in the example you would interpolate your spreadsheet data.

Passing the data that you read from the spreadsheet in as additional parameters rather than reading it in every time the ODE solver calls your ODE function will save time (potentially a lot of time if reading the data takes a while or the ODE solver needs to call your ODE function many times.)

[Bhanu Pratap Akherya]

I did what you said but my function is now returning 0 and NaN values can you tell what could be possibly wrong? sorry, but I am just a beginner.

[Steven Lord]

If you show us the code you've written we may be able to offer suggestions on how to fix and/or improve it.

[Bhanu Pratap Akherya]

Here it is;

main code:

dt=0.001;

T=300;

y0=initial conditions;

[tsol ysol]=ode45('beam_function',[1:dt:T],y0);

plot(tsol,ysol);

function code:

function [dy]=beam_function(t,y)

function F = Force(t)

t_array = xlsread('l&d.xlsx','Q1:Q300'); % This is t array from xls file

f_array = xlsread('l&d.xlsx','O1:O300'); % This is F array from xls file

F = interp1(t_array,f_array,t);

end

dy=[y(1:u-1);

inv(M)*(Force(t)-K*y-C*y)]

end

[Steven Lord]

So this is your main script:

dt=0.001;
T=300;
y0=initial conditions; % I'm assuming this is just abbreviated
[tsol ysol]=ode45('beam_function',[1:dt:T],y0);
plot(tsol,ysol);

My only comment here is that you should pass the first input as a function handle rather than the name of a function. While the syntax you're using does still work for backwards compatibility, the recommended approach for probably close to two decades has been the function handle approach. This also has one more benefit that I'll talk about in a moment.

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

Now for your beam_function function:

function [dy]=beam_function(t,y)
function F = Force(t)
t_array = xlsread('l&d.xlsx','Q1:Q300'); % This is t array from xls file
f_array = xlsread('l&d.xlsx','O1:O300'); % This is F array from xls file
F = interp1(t_array,f_array,t);
end
dy=[y(1:u-1);
    inv(M)*(Force(t)-K*y-C*y)]
end

Every time beam_function calls the Force function nested inside it, it will call xlsread twice to get the same two pieces of data from the spreadsheet. This is inefficient. The anonymous function I passed into ode45 above accepts t_array and f_array as additional parameters, so you can read the data once before calling ode45 and use it over and over. So the updated script and function are:

dt=0.001;
T=300;
y0=initial conditions; % I'm assuming this is just abbreviated
% Read the data from the spreadsheet once
t_array = xlsread('l&d.xlsx','Q1:Q300'); % This is t array from xls file
f_array = xlsread('l&d.xlsx','O1:O300'); % 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);

function [dy]=beam_function(t,y, t_array, f_array)
% Use the data from the spreadsheet over and over
F = interp1(t_array,f_array,t);

% I assume you omitted the code to create variables M, K, and C
% in order to shorten the example
dy=[y(1:u-1);
    inv(M)*(F-K*y-C*y)]
end

[Bhanu Pratap Akherya]

i did exactly as you told but dy is still returning the values as 0 & NaN, my guess is due to the warning shown that was: "Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. " can you tell what is this error and how can I resolve it?

[Steven Lord]

I think we've gotten to the point where debugging that error will be difficult or impossible without your full code (nowhere in this discussion have you shown the code that defines M, K, or C for example) and sample data (synthetic data that "looks like" the real data if you cannot provide the real data.)

[Bhanu Pratap Akherya]

Here is the complete code:

Function code:

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

F = interp1(t_array,f_array,t);

dy=[y(1:u-1);

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

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 = xlsread('l&d.xlsx','Q1:Q300'); % This is t array from xls file

f_array = xlsread('l&d.xlsx','O1:O300'); % This is F array from xls file

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

plot(tsol,ysol(:,Ne))

[Bhanu Pratap Akherya]

in the excel sheet i have bunch of numbers that have been taken from a cfd software.