How to simulate using ode23 with multidimentional second order differential equations

Jesse Crotts (view profile)

on 8 Dec 2018 at 20:20
Latest activity Edited by madhan ravi

on 10 Dec 2018 at 2:43

!!!At the very bottom, there is a much simpler version of this question!!!
I am trying to model a beam using FEM in Matlab. I understand how to use ode23 for single dimentional problems however, my problem is 10 dimentional and second order. It looks like [M]*(xdotdot)+[K]*x=Q where M and K are 10x10 matrices. I think I have my code close to being right but I'm not quite sure what is wrong. When I debug the code it stops when It gets to the ode function at the bottom. Note that it is calling function sys.m
Here is my code:
clc; clear
format shortE
t=1:.01:10; % time
% -------------------------------------------------------
%{
% User input
l = 3; % length in inches
EI1 = 5*10^6; % EI
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = 200; % spring constant 1
k2 = k1; % spring constant 2
m = 268.2; % mass per unit length
po = 100; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = 75; % point load
%}
% -------------------------------------------------------
% -------------------------------------------------------
% User input
% Validation code
l = 3; % length in inches
EI1 = 5*10^6; % EI
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = 0; % spring constant 1
k2 = k1; % spring constant 2
m = 268.2; % mass per unit length
po = 100; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = [75]'; % point load
zeta = .5;
% -------------------------------------------------------
% -------------------------------------------------------
% Phi function handles
phi1 = @(zeta) 1-3*zeta^2+2*zeta^3;
phi2 = @(zeta) l*zeta-2*l*zeta^2+l*zeta^3;
phi3 = @(zeta) 3*zeta^2-2*zeta^3;
phi4 = @(zeta) -l*zeta^2+l*zeta^3;
Phi = @(zeta)[...
phi1(zeta)*phi1(zeta) phi1(zeta)*phi2(zeta) phi1(zeta)*phi3(zeta) phi1(zeta)*phi4(zeta);...
phi2(zeta)*phi1(zeta) phi2(zeta)*phi2(zeta) phi2(zeta)*phi3(zeta) phi2(zeta)*phi4(zeta);...
phi3(zeta)*phi1(zeta) phi3(zeta)*phi2(zeta) phi3(zeta)*phi3(zeta) phi3(zeta)*phi4(zeta);...
phi4(zeta)*phi1(zeta) phi4(zeta)*phi2(zeta) phi4(zeta)*phi3(zeta) phi4(zeta)*phi4(zeta) ];
%PPhi = Phi(zeta)
%}
% This function needs to be moved to where it is needed
% -------------------------------------------------------
% -------------------------------------------------------
% Displacement function handles
Qpo = @(po,l) [po*l/2 po*l^2/12 po*l/2 -po*l^2/12]';
QP = @(zeta,P) [-P*phi1(zeta) -P*phi2(zeta) -P*phi3(zeta) -P*phi4(zeta)]';
Qspring3 = @(k1,zeta) k1*Phi(zeta);
Qspring4 = @(k1,zeta) k1*Phi(zeta);
% This function needs to be moved to where it is needed
% -------------------------------------------------------
% -------------------------------------------------------
% Calculation simplification
a=zeros(4,6);
aa=zeros(4);
aaa=zeros(4,2);
I4=eye(4);
A1=[I4 a];
A2=[aaa I4 aa];
A3=[aa I4 aaa];
A4=[a I4];
% -------------------------------------------------------
% -------------------------------------------------------
% p
% P
% -------------------------------------------------------
% -------------------------------------------------------
% Mass matrix, Stiffness matrix
M1 = (m*l/420)*[...
156 22*l 54 -13*l;
22*l 4*l^2 13*l -3*l^2
54 13*l 156 -22*l
-13*l -3*l^2 -22*l 4*l^2];
M2 = M1; % note, this can be different
M3 = M1; % note, this can be different
M4 = M1; % note, this can be different
K1 = (EI1/l^3)*[...
12 6*l -12 6*l
6*l 4*l^2 -6*l 2*l^2
-12 -6*l 12 -6*l
6*l 2*l^2 -6*l 4*l^2];
K2 = (EI2/l^3)*[...
12 6*l -12 6*l
6*l 4*l^2 -6*l 2*l^2
-12 -6*l 12 -6*l
6*l 2*l^2 -6*l 4*l^2];
K3 = (EI3/l^3)*[...
12 6*l -12 6*l
6*l 4*l^2 -6*l 2*l^2
-12 -6*l 12 -6*l
6*l 2*l^2 -6*l 4*l^2];
K4 = (EI4/l^3)*[...
12 6*l -12 6*l
6*l 4*l^2 -6*l 2*l^2
-12 -6*l 12 -6*l
6*l 2*l^2 -6*l 4*l^2];
% -------------------------------------------------------
% calculations
K3 = K3 + Qspring3(k1,l/d2); % Khat
K4 = K4 + Qspring4(k1,l/d3); % Khat
Qpo4 = Qpo(po,l);
QP4 = QP(l/d1,P);
syms Q11 Q21 Q31 Q41
syms Q12 Q22 Q32 Q42
syms Q13 Q23 Q33 Q43
syms Q14 Q24 Q34 Q44
Q1 = [Q11 Q21 Q31 Q41]';
Q2 = [-Q31 -Q41 Q32 Q42]';
Q3 = [-Q32 -Q42 Q33 Q43]';
Q4 = [-Q33 -Q43 0 0]'+QP4+Qpo4;
Q11 = A1'*Q1;
Q22 = A2'*Q2;
Q33 = A3'*Q3;
Q44 = A4'*Q4;
Q = Q11+Q22+Q33+Q44;
Q(1) = 0; Q(2) = 0;
% -------------------------------------------------------
% -------------------------------------------------------
% global calculations
M = A1'*M1*A1+A2'*M2*A2+A3'*M3*A3+A4'*M4*A4;
K = A1'*K1*A1+A2'*K2*A2+A3'*K3*A3+A4'*K4*A4;
Minv = inv(M);
% -------------------------------------------------------
%qstatic = inv(K)*Q
time = (0:.001:22.5)';
x0 = [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]; %[position,velocity]
[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt]=ode23(@(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt) sys(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt,Minv,K), time, x0);
Then my sys.m looks like:
function [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt] = sys(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt,Minv,K)
f = zeros(20,1);
f(1) = x(2);
f(2) = x(3);
f(3) = x(4);
f(4) = x(5);
f(5) = x(6);
f(6) = x(7);
f(7) = x(8);
f(8) = x(9);
f(9) = x(10);
f(10) = x(11);
delta = -Minv*K*[x(12) x(13) x(14) x(15) x(16) x(17) x(18) x(19) x(20) x(21)]'%+Minv*Q;
f(11) = delta(1)
f(12) = delta(2)
f(13) = delta(3)
f(14) = delta(4)
f(15) = delta(5)
f(16) = delta(6)
f(17) = delta(7)
f(18) = delta(8)
f(19) = delta(9)
f(20) = delta(10)
end
Here is the easier version of the code:
clc; clear
figure
time = (0:.001:22.5)';
M = [1 2; 3 4];
K = [5 6; 7 8];
x0 = [0;0;0;0]; %[position,velocity]
[t1,x1,x2,x3,x4]=ode23(@(t1,x) trick(t1,x,M,K), time, x0);
plot(t1,x(:,1));
xlabel('Time(sec)');
ylabel('Displacement');
title('Stepped Response(Underdamped)');
Here is the code that it calls:
%2x2 sys
function f = trick(t,x,M,K)
f = zeros(4,1);
f(1) = x(2);
f(2) = x(3);
MinvnegK = -inv(M)*K;
delta = MinvnegK*[x(4) x(5)]';
f(3) = delta(1);
f(4) = delta(2)

1 Comment

on 9 Dec 2018 at 8:04
Sorry had to close the other question because this question is the same as the other one .

on 9 Dec 2018 at 9:00

clc; clear
figure
time = (0:.001:22.5)';
M = [1 2; 3 4];
K = [5 6; 7 8];
x0 = [1;0;0;0]; %[position,velocity]
[t1,x]=ode23(@(t,x) trick(t,x,M,K), time, x0);
plot(t1,x(:,1));
xlabel('Time(sec)');
ylabel('Displacement');
title('Stepped Response(Underdamped)');
%2x2 sys
function f = trick(t,x,M,K)
f = zeros(4,1);
f(1) = x(1);
f(2) = x(2);
MinvnegK = -inv(M)*K;
delta = MinvnegK*[x(3) x(4)]';
f(3) = delta(1);
f(4) = delta(2);
end

Jesse Crotts

Jesse Crotts (view profile)

on 10 Dec 2018 at 2:13
Thanks for the answer. Do you know of a faster way to write the values at the end of the function. I know its easy for a 2x2 matrix but if you had a 100x100 matrix what would be the best way to do it:
delta = MinvnegK*[x(3) x(4)]';
f(3) = delta(1);
f(4) = delta(2);