ode45 with matrix initial conditions

Question

I am a novice user of MATLAB. I have code that is trying to find the time history of a state space model. There are four first order ODEs that I want to solve simultaneously using ode45. The essence of the equations to be solved is as follows:

x1_dot = x2

x2_dot = -[M] * [K] * x1 - [M] * [C] * x2 + constant*[M] * [P3] * x3 + constant*[M] * [P4] * x4

x3_dot = x2 - constant*x3

x4_dot = x2 - constant*x4

Where [M], [K], [C], [P3], and [P4] are 3x3 matrices; x1, x2, x3, x4 are all 3x1 vectors; and x1_dot etc. represent the time derivatives (which are 3x1 vectors). I have initial conditions only for x1.

The MATLAB code I've written is below. This code is within my overall program. I am not calling a separate function because I do not know how to pass all of the matrices/vectors into ode45 through a function. I am getting the error: "Index exceeds matrix dimensions."

tspan = 0:1:20;

initial = [0 0.03491 0];

f = @(t,x) [x(2);
            -inv(M_Dbl_Bar_Matrix)*K_Dbl_Bar_Matrix*x(1) - inv(M_Dbl_Bar_Matrix)*C_Dbl_Bar_Matrix*x(2) + (0.5*rho*U^2)*inv(M_Dbl_Bar_Matrix)*P3_Matrix*x(3) + (0.5*rho*U^2)*inv(M_Dbl_Bar_Matrix)*P4_Matrix*x(4);
            x(2) - Beta_1*x(3);
            x(2) - Beta_2*x(4)];

[t,xp] = ode45(f,tspan,initial);

Questions:

How do I address the x(1), x(2), x(3), and x(4) in ode45 being 3x1 vectors?
How do I apply initial conditions for this system of equations? For example, do I use a vector format such as: initial = [0 0.03491 0; 0 0 0; 0 0 0; 0 0 0]?
Am I correctly using / written the function (f) and ode45?

Rody Oldenhuis Rody Oldenhuis · Accepted Answer · 2017-01-05T09:45:52

Just integrate the 12×1 system as if they are 12 coupled ODEs.

A couple of other observations:

Avoid inv() wherever possible -- it's slow and inaccurate. Use mldivide (or mrdivide) instead. Moreover, you recompute it no less than 4 times at each evaluation of f, while it's basically constant!
You seem to want output at every second. ode45 is a variable step integrator, meaning that it automatically adjusts the step size based on estimates of the error made at each step. Requesting outputs at times that deviate from the times chosen by ode45 (called "dense output") is easy to do, but it'll cost you extra function evaluations. Usually that's not really needed, and you can get away with using the more efficient tspan = [0 20] instead. It all depends on your specific needs.

Now, here's what I came up with:

% Time interval of interest
tspan = [0 20];

% Initial values
x1_0 = [0 0.03491 0];
x2_0 = [0       0 0];
x3_0 = [0       0 0];
x4_0 = [0       0 0];
x0   = [x1_0 x2_0 x3_0 x4_0].';

% Pre-compute a few constants
Fd  = 0.5*rho*U^2;
P3f = Fd*M_Dbl_Bar_Matrix\P3_Matrix;
P4f = Fd*M_Dbl_Bar_Matrix\P4_Matrix;
P1f = -M_Dbl_Bar_Matrix\K_Dbl_Bar_Matrix;
P2f = -M_Dbl_Bar_Matrix\C_Dbl_Bar_Matrix;

% The derivative (12×1, but constructed as 4·3×1)

one = 1:3;   three = 7:9;   % well-named index vectors to
two = 4:6;   four  = 10:12; % make our lives a bit easier

f = @(t,x) [x(two)
            P1f*x(one) + P2f*x(two) + P3f*x(three) + P4f*x(four)
            x(two) - Beta_1*x(three)
            x(two) - Beta_2*x(three)];

Now integrate this system
[t, xR] = ode45(f, tspan, x0);

% extract results in the same kind of blocks:
x1 = xR(:, one);
x2 = xR(:, two);
x3 = xR(:, three);
x4 = xR(:, four);

% ... process the results in whatever way you see fit

ode45 with matrix initial conditions

2 Answers

EDIT