how to use ode45 to solve motion equation with Matrix

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')

[alize beemiel]

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

Answer 1

[Alan Stevens]

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
   

[alize beemiel]

oh !!!

now its clear ..

thank you very much sir