# how to solve matrix differential equation

9 views (last 30 days)
riccardo carrieri on 19 Aug 2019
Edited: David Wilson on 20 Aug 2019
hi,
I' m trying to solve this differential equation for education purposes but I'm having difficulty to solve it.
does someone know a simple way to solve it?
? ∗ ?̈ = ?? − (?? + ??) ∗ ?
I attached a pdf where you can find all the information about the program.
thank you very much.

David Wilson on 20 Aug 2019
Edited: David Wilson on 20 Aug 2019
Here's my attempt:
M = [
97227/21875, 0, 82269/109375, 0, 67311/43750, 0, -97227/218750, 0; ...
0, 97227/21875, 0, 82269/109375, 0, 67311/43750, 0, -97227/218750;
82269/109375, 0, 89748/546875, 0, 97227/218750, 0, -67311/546875, 0;
0, 82269/109375, 0, 89748/546875, 0, 97227/218750, 0, -67311/546875;
67311/43750, 0, 97227/218750, 0, 97227/21875, 0, -82269/109375, 0;
0, 67311/43750, 0, 97227/218750, 0, 97227/21875, 0, -82269/109375;
-97227/218750, 0, -67311/546875, 0, -82269/109375, 0, 89748/546875, 0;
0, -97227/218750, 0, -67311/546875, 0, -82269/109375, 0, 89748/546875]
Qv = [ 0; -20651534751534355/422212465065984; 0; -4130306950306871/422212465065984; 0;
-20651534751534355/422212465065984; 0; 4130306950306871/422212465065984];
I couldn't be bothered using your matrices, so I just used some random, symmetric ones. You'll need to type in yours ...
Kt = randn(8,8); Kt = Kt'*Kt; % force symmetric
Kl = randn(8,8); Kl = Kl'*Kl; % force symmetric
L = 5; % ???
e0=[0; 0; 1; 0; L; 0; 0; 1];
Now I'm going to
(1) Convert the 2nd order ODE to Cauchy form, and
(2) Generate the appropriate Mass matrix, (although you can easily invert this since the condition number is reasonable.)
Once all set up, just throw to ode45 to solve it. Once done, extract the first 8 elements & plot them.
%-------------------------------------------------
K = (Kt+Kl);
Ma = eye(16); % augmented mass matrix.
Ma(9:end, 9:end) = M; % Bottom right corner is "M"
z0 = [e0; zeros(8,1)];
odefun = @(t,z) [z(9:16,1); ...
Qv - K*z(1:8,1)];
tspan = [0,5]; % who knows what this should be ???
optns = odeset('Mass',Ma); % Include Mass matrix & keep on LHS.
[t,z] = ode45(odefun, tspan, z0, optns );
plot(t,z(:,[1:8]));
##### 2 CommentsShowHide 1 older comment
David Wilson on 20 Aug 2019
First up, you shouldn't use the inv(Ma) as well as the "Mass" matrix. Keep the mass matrix on the LHS, and drop the inv(Ma) part.
If you are going to use the inv(Ma), then drop the mass matrix part, however this is not generally recommended nor efficient.
You error seems to suggest that one (or more) ofthe arguments to ode45is the wrong type. Check each of the arguments.