plotting matlab state space

plotting matlab state space

I am trying to simulate a very simple dynamical system.

mx''+bx'+kx=u(t)
(mx''+cx'+kx=m*9.81)

I want to draw the displacement (position x) and velocity (v) of the system over a time from 0 to 10s with 3 different time stepsizes. i am supposed to work with an A-matrix approach;

ydot = [(-c/m)(-k/m);1 0] [x xdot] +[9.81;0]

i have no idea where to start due to the fact im new to matlab. Here is what i have

clc;
close all;
clear all;
 
m=7;
c=10;
k=100;
dt=5
t=0:dt:10
A=[(-c/m) (-k/m);1 0 ];
B=[9.81 ; 0];

Any suggestion for completing the code?

Best regards,

Marco

[darova]

I know how to solve this using ode45 (you can get x and

mx''+bx'+kx=u(t)

But don't know what is this for

ydot = [(-c/m)(-k/m);1 0] [x xdot] +[9.81;0]

[Marco Boutrus]

I need to write a script for the rewritten form in state space.

Form: ydot = Ay + f where A is the system matrix and vector f represents the inhomogeneous term. Which is given:

ydot = [(-c/m)(-k/m);1 0] [x xdot] +[9.81;0]

And plot the results for the position and the velocity of the mass for t=0 to 10s with three different time steps

[darova]

as i understood you correctly

[Marco Boutrus]

let me quote the assignment:

4. Solving numerically:

o Rewrite the equation in state-space form, using the state vector ? = (?). The ?̇

end result should look like this:

?̇ = A? + ?

where A is called the system matrix and vector f represents the inhomogeneous term. Find the expressions of the matrix A and vector f. Note that not all terms need to be present, i.e. some matrix or vector elements may be zero.

o Write a Matlab script that solves the equations of motion numerically for the parameter values given to you (group dependent).

N.B.: Use the A-matrix approach as discussed. You will need the script later on for a different system, having a different A-matrix. This way you can re-use your Matlab code.

o Determine the solution of the system above and plot the results for the position and velocity of the mass for t = 0 to 10s using the three different time steps given in the parameter set.

5. Verification/validation

o Check the numerical solution for the different time steps with the analytical

solution. What do you see? Which of the solutions is/are the correct one(s)?

To gain some more insight into the behaviour of the system, we can have a look at the eigenvalues of the A-matrix. This yields valuable information about the behaviour of the system, even before the solution is determined.

o Determine the eigenvalues of the A-matrix and discuss these in comparison to the solutions you have plotted. Check the real and imaginary parts of the eigenvalues and compare those with the results you plotted.

▪ What does the real part of the eigenvalues do to your solution?

▪ What does the imaginary part of the eigenvalues do to the solution? ▪ Measure the amplitude and period of the solution and discuss these

values in relation to the eigenvalues.

o What can you say about the contribution of the inhomogeneous term in the

solution?

Answer 1

[Ameer Hamza]

Following use ode45 to solve your ODE

dt=5;
t= [0 10];
ic = [0;0]; % initial condition

[t,y] = ode45(@myOdeFcn, t, ic);
plot(t,y);
legend({'Velocity', 'Position'});

function dydt = myOdeFcn(t, y)
    m=7;
    c=10;
    k=100;
    A=[(-c/m) (-k/m);1 0 ];
    B=[9.81; 0];
    dydt = A*y(:)+B;
end

[Marco Boutrus]

Thank u so much for your answer, could u care to elaborate the meaning of the 0defcn? also when im trying to change the stepsize, the graph does not change

[Ameer Hamza]

Marco, myOdeFcn is just the implementation of the ODE you wrote in your question. You can see it has the same equation as yours.

About step-size, I didn't understand your point. You you just what to make graph at 3 points t=[0 5 10]. That will appear a bit strange but if that is what you want then try

dt=5;
t= 0:dt:10; % <---- this line is changed
ic = [0;0]; % initial condition

[t,y] = ode45(@myOdeFcn, t, ic);
plot(t,y);
legend({'Velocity', 'Position'});

function dydt = myOdeFcn(t, y)
    m=7;
    c=10;
    k=100;
    A=[(-c/m) (-k/m);1 0 ];
    B=[9.81; 0];
    dydt = A*y(:)+B;
end

[Marco Boutrus]

ahaa i get it now, however its not exactly what im looking for. i did post a small text above with the details of what i need exactly! if u can and got time, it would help me loads if ud take a look at it

[Ameer Hamza]

I saw that comment. I think that your professor what to write your own ode solver instead of using ode45. Probably he has taught a solver in his lectures, which he wants you to use. Without the formula of the solver given by the professor, it is difficult to suggest a solution.

[Marco Boutrus]

Yes that is true, however due to the corona virus we didnt get the lecture where the professor wouldve taught us that, which is the reason why im here :( i dont think he is only taking one way of an answer. if u can explain what is happening usually thats fine

[Marco Boutrus]

i have calculated the formula to this form ?̇ = A? + ? by hand on paper -> ?̇ = [(-c/m)(-k/m);1 0] [x xdot] +[9.81;0], but the writing a script part is the hard part. i think myself, it has something to do with either the Euler or Heun scheme

[Ameer Hamza]

Yes, you are correct. You will need to write the code for something like that. Since this is a homework question, you first need to write the code for the algorithm yourself, and we can help you if you face some issues, specifically related to MATLAB.

[Marco Boutrus]

Yes thank you! could u give me a starting hint? hahah

[Ameer Hamza]

See here: https://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b#answer_107643

Few of the MATLAB's old ode solvers were based on Euler and Runge-Kutta. MATLAB have still made their code available. Download and study them. They will give you the general structure of the numerical solver. Adapt it according to your function.

[Marco Boutrus]

When following the heun method as given by the professor, this is what i get. however its showing the graph reversed(starting left instead of right) and only the Displacement is shown. where is the fault?

clear all
clc
close all
 
m=7
c=10
k=100
% Define the system matrix of the differential equation
A=[(-c/m) (-k/m);1 0 ];
B=[9.81;0]
%dydt= A*y+B
 
dt = 0.5;   % Time step [s]
t = 0:dt:10; % Time range [s]
 
y = zeros(2, length(t)); % Pre-initialise the solution vector y = [x;x']
 
x0     = 1;       % Initial displacement [m]
v0     = 0;       % Initial velocity x' [m/s]
y(:,1) = [x0;v0]; % Initial conditions, first column of solution vector
 
% Prepare the plotting window
figure(1)
hold on
 
% Numerical integration using Heun scheme
for i=1:length(t)-1
    y_star = y(:,i) + A*y(:,i)*dt;                    % Predictor
    y(:,i+1) = y(:,i) + 0.5*(A*y(:,i) + A*y_star)*dt; % Corrector
end
plot(t,y(1,:),'b--') 
 
% Exact solution to the differential equation
% The eigenvalues are complex conjugates, so a damped oscillation
t_e = 0:0.5:10; % Define new time vector (fine grid for smooth lines)
E = eig(A);      % Eigenvalues of the differential equation
% Determine real and imaginary parts of the two eigenvalues
alfa = real(E(1));
beta = abs(imag(E(1)));
% Initial conditions solve for the coefficients of the solution
C1 = x0;
C2 = (v0-alfa*C1)/beta;
% Exact solution
x_e = exp(alfa*t_e).*(C1*cos(beta*t_e) + C2*sin(beta*t_e));
plot(t_e,x_e,'r-')
 
grid on
xlabel('Time [s]')
ylabel('Displacement [m]')
legend('Heun','Exact')

[Ameer Hamza]

This happens because you are using too large value of dt=0.5. The solver diverges. Use a small value say dt=0.1 and see the result.

[Ameer Hamza]

Also, if you want to show both displacement and velocity then change the first plot line like this

plot(t',y','b--') 

[Marco Boutrus]

but i am supposed to see the graph with dt also as 0.5? i have changed those small edits, but now the v graph doesnt use v0 as ic. its very confusing hahah

also atm im not using the B matrix? and im only using the A matrix to solve the eigenvalues, is that correct?

And if i change x0 to 0, the whole activity is flatlined

[Ameer Hamza]

I guess the 0.5 second interval was required to demonstrate that the solver can diverge if dt is large. Also, when I plot with

plot(t',y','b--') 

both x and v follow the initial conditions.

I am not sure about the B matrix, but I guess you can replace A*y with A*y+B to solve the omplete equation.

Since you are just using A*y so this is the natural response of the system. Therefore, If the initial position of the system is x=0, then you don't see any oscillations.

[Marco Boutrus]

ah yes now all of sudden it does show me both starting from the initial conditions. however i think im changing the wrong part to ay+b, which one needs to be changed? also is the eigenvalue part correct?

[Ameer Hamza]

The eigenvalue part is not useful for the numerical solution. The first part should be enough for your assignment. Matrix B can be included like this

for i=1:length(t)-1
    y_star = y(:,i) + (A*y(:,i)+B)*dt;                    % Predictor
    y(:,i+1) = y(:,i) + 0.5*(A*y(:,i) + B + A*y_star + B)*dt; % Corrector
end

Also note that due to the way matrix A is created, the first variable is velocity, and the second variable is the position. Take care of this when interpreting the variables

x0     = 0;        % Initial displacement [m]
v0     = 1;        % Initial velocity x' [m/s]
y(:,1) = [v0; x0]; % Initial conditions, first column of solution vector

[Marco Boutrus]

Yes ive changed those small things. but somehow when i change dt from 0.05 and 0.1 no changes are viewable in the graph. the graph stays the same. also when i change the dt to 0.5, for some reason the graph starts on the right side instead of the left.

also the exact solution is just a line through the middle, im guessing because of the fact that part B is not incorporated:S how can i fix that?

Its so confusing for a starter

[Ameer Hamza]

Marco, when dt=0.5, the graph does not start from the right. The graph is actually divergeing. It means that it will continue to increase to infinity. You can increase the end time t=[0 100] and see the y-axis. It will show that the solution is diverging.

[Marco Boutrus]

but looking at the formula, how can it not start at 0?

And how come when dt is changed between 0.05 and 0.1, the graph stays exactly the same?

[Ameer Hamza]

I am not sure what you meant by the graph not starting from zero.

The graph is not exactly same for 0.05 and 0.1. But the difference is too minor to be observed with eyes. You can check the number of elements of t and y to see the difference for both cases.

[Marco Boutrus]

ahh yes it was indeed to minor to see with the naked eye.

But i still cant put my finger on the dt=0.5 graph. it goes to an insane displacement and velocity. how can dt change the whole graph in this way?

[Ameer Hamza]

This is nothing strange, and it is a well-known fact. All numerical method suffers from this to some extent. Note that Euler is one of the basic numerical methods, so it will be affected the most. See the numerical stability section here: https://en.wikipedia.org/wiki/Euler_method#Numerical_stability