changing values of RHS with each time step in ODE

Question

Bhanu Pratap Akherya 2021-8-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1436182-changing-values-of-rhs-with-each-time-step-in-ode

评论： Wan Ji 2021-8-24

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?

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

请先登录，再进行评论。

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

Answer 1

Wan Ji 2021-8-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1436182-changing-values-of-rhs-with-each-time-step-in-ode#answer_769702

在 MATLAB Online 中打开

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

Wan Ji 2021-8-20

在 MATLAB Online 中打开

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

编辑：Wan Ji 2021-8-20

在 MATLAB Online 中打开

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

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

Wan Ji 2021-8-24

在 MATLAB Online 中打开

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

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1436182-changing-values-of-rhs-with-each-time-step-in-ode#answer_769767

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

Steven Lord 2021-8-18

在 MATLAB Online 中打开

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

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

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

请先登录，再进行评论。

changing values of RHS with each time step in ODE

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

回答（2 个）

9 个评论
显示 7更早的评论隐藏 7更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

changing values of RHS with each time step in ODE

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

回答（2 个）

9 个评论 显示 7更早的评论隐藏 7更早的评论

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论