how to use ode45 to solve motion equation with Matrix

2 次查看（过去 30 天）
alize beemiel
alize beemiel 2023-4-1
评论： alize beemiel 2023-4-1
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 个评论
alize beemiel
alize beemiel 2023-4-1
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

请先登录，再进行评论。

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

采纳的回答

Alan Stevens
Alan Stevens 2023-4-1
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by