how to use ode45 to solve motion equation with Matrix

1 Apr 2023
1 Answer
Answer Accepted
Updated 1 Apr 2023
8 Views (30 days)

0 votes

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
Show -1 older comments Hide -1 older comments

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.

Sign in to answer this question.

 Accepted Answer

Alan Stevens
Alan Stevens on 1 Apr 2023
Ran in:

1 vote

Ran in:
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 Ordinary Differential Equations in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!