How to simulate using ode23 with multidimentional second order differential equations
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Jesse Crotts
am 8 Dez. 2018
Bearbeitet: Bob
am 13 Feb. 2023
!!!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];
% -------------------------------------------------------
% -------------------------------------------------------
% load transformations
% 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 Kommentar
madhan ravi
am 9 Dez. 2018
Bearbeitet: madhan ravi
am 9 Dez. 2018
Sorry had to close the other question because this question is the same as the other one .
Akzeptierte Antwort
madhan ravi
am 9 Dez. 2018
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
![Screen Shot 2018-12-09 at 2.25.48 PM.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/198606/Screen%20Shot%202018-12-09%20at%202.25.48%20PM.png)
2 Kommentare
madhan ravi
am 10 Dez. 2018
Bearbeitet: madhan ravi
am 10 Dez. 2018
delta(:,100) just use indexing to store the values
Weitere Antworten (1)
Bob
am 13 Feb. 2023
Bearbeitet: Bob
am 13 Feb. 2023
It is late :( but might be useful :)
clear; clc; % clear the desk
t = [ 0 , 10 ] ; % time interval, use the default interior points
x = [1 0, 0 0]'; % starting point [position, velocity]
B = inv([1 2; 3 4]) * [5 6; 7 8] ; % define matrix B = -inv(M) * K
s = ode23 (@(tau,z) [z(1:2); B*z(3:4)], t, x); % solve with ode23()
plot(s.x, s.y(1,:));% check the shape of the solution and its boundaries
disp(sprintf("\tt in [ %d, %d ] \ts in [ %g, %g ]",t,min(s.y(1,:)), max(s.y(1,:))));
0 Kommentare
Siehe auch
Kategorien
Mehr zu Matrices and Arrays finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!