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:

1. How do I address the x(1), x(2), x(3), and x(4) in ode45 being 3x1 vectors?

2. 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]?

3. Am I correctly using / written the function (f) and ode45?

Answer 1

1.) What do you mean by "being 3x1 vectors"? x(i) has to be a single variable in order for this to work, so size(x)=1x4. Having x(i)=(x,y,z) does not really make any sense in the context of ODE's.

2.) Your variable vector x should have length 4, so your initial conditions should reflect this.

initial = [0 0.03491 0 0];

Works fine for me.

3.) Yes, I think so.

Could you maybe provide more information on what you are trying to do?

EDIT

Ok, I think I understand now what you are trying to do.

You have 4 vectors for instance: x(x,y,z), x'(x,y,z), p(x,y,z), q(x,y,z) and a large matrix 4x4 matrix consisting of 3x3 matrices. So I guess

(forgive the hosted images, since im a stackoverflow noob I'm not allowed to post html or images directly, or multiple links)

Reference for Math 1,2,3: https://postimg.org/gallery/1qh2ywiqq/

Math 1 from link

with

Math 2 from link

your state space equation.

So to solve this you have to set up a 12 dimensional ode45 problem, instead of the 4 dimension problem you have now, since each vector has 3 components. What you have to do is change the 4x4 matrix into a 12x12 matrix by specifying each of the entries explicitly. You also have to give f = @(t,x) a 1x12 vector:

Math 3 from link

Then set initial to a 1x12 vector of initial conditions (one for each dimension in our vector).

With the matrix form and initial conditions we have now, we can use this source: https://nl.mathworks.com/help/symbolic/solve-a-system-of-differential-equations.html#buxuujb

which tells us how to set this up properly.

I currently don't have the time to give you the code, but I hope I understood correctly what you are trying to do and this will be useful.

Answer 2

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

A couple of other observations:

1. 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!
2. 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