How to use Runge Kutta method to solve state space equation

I am trying to use the Runge Kutta method to solve

xdot = Ax + Bu

Matrices A and B are constants and the input vector u is a known function of time. I am looking to get values for x. I have seen lots of posts about runge kutta but am very confused as to how to get it to work in matlab.

[James Tursa]

What have you done so far? What specific problems are you having with your code?

[Eilidh Graham]

% Solve the equations of motion using a Runge-Kutta routine
[t,x2]=ode45(@LinearEoM,t,uc);

save (fname, 't', 'x2')


function xdot2 = LinearEoM(t, x2)
global Ac Bc uc 
    xdot2 = Ac*x2 + Bc*uc;
end

Ac is a 4x4 matrix, Bc 4x12, uc is a timeseries with 12 rows. I am looking to get x2 as a timeseries, it will have 4 rows. I am just not sure if i am doing this right at all

Answer 1

[Walter Roberson]

(4 by 4) * (4 by 1) --> 4 by 1

(4 by 12) * (12 by 1) --> 4 by 1

4 by 1 + 4 by 1 --> 4 x 1

So that calculation sounds plausible -- but only if your uc timeseries has only 1 column.

However... you are passing in uc as the initial conditions, as well as using uc independently inside the code. If uc has 12 rows and it also forms the initial conditions, then the first x2 value passed into the code would have 12 rows, and you would be ding (4 by 4) * (12 x something) which is going to fail.

When you are looking to get an output that has 4 rows, then your initial conditions should have 4 rows.

Also we recommend against using global. See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html . Oh and remember that Bc*uc is independent of your x2 value, so there is no point re-calculating Bc*uc each iteration: might as well calculate it once outside the function and pass the combination in.

[Eilidh Graham]

Many thanks, I think I have got a bit further now

%initial conditions for the integration (the trim state)
x2_0 = [0;0;theta_new;phi_new];
% Solve the equations of motion using a Runge-Kutta routine
[time,x2]=ode45(@LinearEoM,[0 t1],x2_0);

then my function is (still need to try not use global)

function xdot2 = LinearEoM(time, x2)
    global Ac Bc uc
    xdot2 = Ac*x2 + Bc*uc;
end

I am still a bit confused about the uc bit - as a timeseries it is a 1x1timeseries but if it is just a matrix it is a 12x# matrix. The number of columns depends on a time chosen. will i need to run it for howver many columns there are. thanks in advance

[James Tursa]

Can you give the exact dimensions of uc? Normally I would have expected a need to interpolate uc to the proper time before multiplying by Bc.

[Eilidh Graham]

For the time value i am using currently, it is a 12x21 matrix. Been playing around with loops try but not sure if that i sthe right approach

[James Tursa]

But if Bc is 4x12 and uc is 12x21, the product Bc*uc would be 4x21 which isn't going to work. Hence my question about interpolation. Are the columns of uc the controls at different times? Are you sure you don't want to interpolate the 12x21 uc matrix to a single 12x1 vector that matches time?

[Eilidh Graham]

Yes, the columns of uc are the control values at different timepoints. How would I go about interpolating the uc matrix to make it into a 12x1 vector? Many thanks

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/interp1.html#btwp6lt-1-v

Notice the part where you can pass in the samples as a 2d array. You can take advantage of that to interpolate on each row of the 12 x 21 separately.

Be sure to use 'spline' or similar interpolation method. Do not use the default linear interpolation. The mathematics of ode*() requires that the calculations be (continuous second derivatives), which spline will give you but linear cannot.

[Eilidh Graham]

Sorry, I am still confused as to how to interpolate my uc matrix. I need x2 to be coming out as a 4x# of timepoints so as i am able to plot it

[Walter Roberson]

You indicate that your time series has 12 rows and 21 columns. Timeseries have one row per time instance so that would describe a time series with 12 time instances and 21 elements for each time. But your code requires that you find 12 values at any given time, which would imply that the columns are time and the rows are the values you need.

To be more correct, timeseries objects can have two different organizations:

1. 2 dimensional data in which number of rows matches number of times, and columns are the data; OR
2. 3 or more dimensional data in which the final dimension must match the number of times, and the previous dimensions are the data for that time

As you have 2 dimensional data that is 12 x 21 then the number of time instances for you data must be 12. But the structure of your problem implies that the number of time instances should be 21 and that you want 12 values at each time.