force input in ode

Question

Bhanu Pratap Akherya 2021-8-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1439834-force-input-in-ode

编辑： darova 2021-8-29

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(u:end)-C*y(1:u-1))]

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?

Note: F here is not time dependent rather it has been taken from results obtained from a simulation software. Then what should be my approach?

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

请先登录，再进行评论。

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

Answer 1

darova 2021-8-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1439834-force-input-in-ode#answer_773639

在 MATLAB Online 中打开

You need to pass F into your ode function. Read more: LINK

for i = 1:length(F)
    [t,y] = ode45(@(t,y)beam_function(t,y,F(i)),...)
end
function dy = beam_funciton(t,y,F)

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

Bhanu Pratap Akherya 2021-8-24

编辑：darova 2021-8-25

在 MATLAB Online 中打开

that helped, i have one more question

here is my main code:

clear all
clc
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=1/(max_freq*20);
dt=0.001;
T=300;
%Displacement initials
y0=zeros(2*(2*(Ne+1))-4,1);
% y0(1:2,:)=[];
% y0(u:u+1,:)=[];
y0(end-1,1)=0.5;
% t_array = a(1,:); % This is t array from xls file
F=xlsread('l&d.xlsx','O1:O300'); % This is F array from xls file
%% ode solver
for i = 1:length(F)
[tsol ysol]=ode15s(@(t, y) beam_function(t, y, F(i), K, M, C, u),[1:dt:T],y0);
end
ysol(:,u:end)=[];%eliminate velocity nodes 5001*12
size(ysol)
 ii=1:2:(u-2); %eliminate angular displacement nodes
  ysol(:,ii+1)=[]; %eliminating angular dis components
  size(ysol)
  
%% plotting
plot(tsol,ysol(:,Ne))

here is my function

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

I want force (F) to be applied on the last node. Can you tell me a way how i can do that?

Bhanu Pratap Akherya 2021-8-25

No, I meant how to apply the current F matrix on only the last column of y in ODE. Right now it is being applied on all the columns of y.

darova 2021-8-29

编辑：darova 2021-8-29

在 MATLAB Online 中打开

Maybe you need another ode function

function [dy]=beam_function(t,y)
n = numel(y);
dy(n/2) = y(end);
dy(end) = inv(M)*(F-K*y(end)-C*y(end));

And solve it separately

请先登录，再进行评论。

force input in ode

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

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

force input in ode

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

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

WeChat

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

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