Clear Filters
Clear Filters

how to use ode45 to solve motion equation with Matrix

8 views (last 30 days)
Hi every one
I want to use ode45 for solving motion equatiom
the equation is a second_oder_ode
% M * (Z)'' + R*( Z)' + K *(Z) = 0
the unknow is Z
My code is
dt=0.1;
t_ode=0:dt:10;
Z0 =[ 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0];
M =[ 0.0210 0.0014 -0.0209 0 0 0 0 0
0.0014 8.1991 -0.0000 -3.9745 -0.0014 0 0 0
-0.0209 -0.0000 0.0420 0.0014 -0.0209 0 0 0
0 -3.9745 0.0014 8.1991 -0.0000 -3.9745 -0.0014 0
0 -0.0014 -0.0209 -0.0000 0.0420 0.0014 -0.0209 0
0 0 0 -3.9745 0.0014 8.1991 -0.0000 -0.0014
0 0 0 -0.0014 -0.0209 -0.0000 0.0420 -0.0209
0 0 0 0 0 -0.0014 -0.0209 0.0210];
K =[ 0.1808 -0.9600 0.0592 0 0 0 0 0
-0.9600 23.3600 0.0000 -11.6800 0.9600 0 0 0
0.0592 0 0.3617 -0.9600 0.0592 0 0 0
0 -11.6800 -0.9600 23.3600 0.0000 -11.6800 0.9600 0
0 0.9600 0.0592 0 0.3617 -0.9600 0.0592 0
0 0 0 -11.6800 -0.9600 23.3600 0.0000 0.9600
0 0 0 0.9600 0.0592 0 0.3617 0.0592
0 0 0 0 0 0.9600 0.0592 0.1808];
R =[ 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 ] ;
odefun = @(t,Z) [Z(2);
-(R*Z(2)+K*Z(1))/M]
[T,Z] = ode45(odefun,t_ode,Z0')
  1 Comment
alize beemiel
alize beemiel on 1 Apr 2023
someone give me this solution
but i dont understaiding
function ode_test
clear all
clc
dt=0.1;
t_ode=0:dt:10;
Z0 =[0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
M =[ 0.0210 0.0014 -0.0209 0 0 0 0 0
0.0014 8.1991 -0.0000 -3.9745 -0.0014 0 0 0
-0.0209 -0.0000 0.0420 0.0014 -0.0209 0 0 0
0 -3.9745 0.0014 8.1991 -0.0000 -3.9745 -0.0014 0
0 -0.0014 -0.0209 -0.0000 0.0420 0.0014 -0.0209 0
0 0 0 -3.9745 0.0014 8.1991 -0.0000 -0.0014
0 0 0 -0.0014 -0.0209 -0.0000 0.0420 -0.0209
0 0 0 0 0 -0.0014 -0.0209 0.0210]
K =[ 0.1808 -0.9600 0.0592 0 0 0 0 0
-0.9600 23.3600 0.0000 -11.6800 0.9600 0 0 0
0.0592 0 0.3617 -0.9600 0.0592 0 0 0
0 -11.6800 -0.9600 23.3600 0.0000 -11.6800 0.9600 0
0 0.9600 0.0592 0 0.3617 -0.9600 0.0592 0
0 0 0 -11.6800 -0.9600 23.3600 0.0000 0.9600
0 0 0 0.9600 0.0592 0 0.3617 0.0592
0 0 0 0 0 0.9600 0.0592 0.1808]
R =[ 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 ]
[T,Z]=ode45( @(t,Z) ode_1Dr(Z,M,R,K) ,t_ode,Z0' )
end
function [dZl] = ode_1Dr (Z,M,R,K)
B=[M zeros(length(M))
zeros(length(M)) M];
C=[R K
-M zeros(length(M))];
dZl=B\(-C*Z);
end

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 1 Apr 2023
Does this help?
%% Data
M =[ 0.0210 0.0014 -0.0209 0 0 0 0 0
0.0014 8.1991 -0.0000 -3.9745 -0.0014 0 0 0
-0.0209 -0.0000 0.0420 0.0014 -0.0209 0 0 0
0 -3.9745 0.0014 8.1991 -0.0000 -3.9745 -0.0014 0
0 -0.0014 -0.0209 -0.0000 0.0420 0.0014 -0.0209 0
0 0 0 -3.9745 0.0014 8.1991 -0.0000 -0.0014
0 0 0 -0.0014 -0.0209 -0.0000 0.0420 -0.0209
0 0 0 0 0 -0.0014 -0.0209 0.0210];
K =[ 0.1808 -0.9600 0.0592 0 0 0 0 0
-0.9600 23.3600 0.0000 -11.6800 0.9600 0 0 0
0.0592 0 0.3617 -0.9600 0.0592 0 0 0
0 -11.6800 -0.9600 23.3600 0.0000 -11.6800 0.9600 0
0 0.9600 0.0592 0 0.3617 -0.9600 0.0592 0
0 0 0 -11.6800 -0.9600 23.3600 0.0000 0.9600
0 0 0 0.9600 0.0592 0 0.3617 0.0592
0 0 0 0 0 0.9600 0.0592 0.1808];
R =[ 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 ] ;
G = [zeros(8) ones(8);
M\(-R) M\(-K)];
%% Initial conditions
Z0 =[ 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0];
% The first 8 values are displacements, the last 8 are velocities.
%% Times
dt = 0.1;
tspan = 0:dt:10;
%% ODE call
[t,Z]=ode45(@(t,Z) odefun(t,Z,G),tspan,Z0);
X = Z(:,1:8);
V = Z(:,9:16);
%% Plots
subplot(2,1,1)
plot(t,X),grid
title('displacements')
xlabel('time'),ylabel('X')
subplot(2,1,2)
plot(t,V),grid
title('velocities')
xlabel('time'),ylabel('V')
%% Function
function dZdt = odefun(~,Z,G)
X = Z(1:8); % current displacements
V = Z(9:16); % current velocities
dZdt = G*[X; V];
end

More Answers (0)

Categories

Find more on Dialog Boxes in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!