Use Euler's method for Mass-Spring System

Hello everybody.

I need to implement Euler's method on a equation based in Mass-Spring System which is:

(m((d^2)x)/(d(t^2)))+(c(dx/dt))+kx=0

Where my x is the displacement (meters), t is the time (seconds), m the mass which is stated as 20kg, my c=10, is the cushioning coefficient and k is the spring value of 20N/m.

So, as my inicial x=1, I need to solve this by Euler with the time interval between 0<=t<=15 seconds.

And them plot the results.

Can you guys help me with the code to solve this?

Thanks so much in advance.

Sander

[James Tursa]

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

Hint: Solve your DE for the highest derivative. E.g., finish this line:

(d^2)x/d(t^2) = ____________

You have a 2nd order system, so you will be integrating two variables at each step, x and dx/dt. Write out derivative equations for each, and then that is what you will use in your Euler step.

[Sander Z]

Is there any example that I can view how to implement a code for this type of equation? I'm kinda lost on where to begin with.

[James Tursa]

You can read here for starters. It even has a MATLAB code example for one variable (but your case will have two variables):

https://en.wikipedia.org/wiki/Euler_method

Answer 1

[James Tursa]

I will get you started. Let y be a 2-element vector containing your states. Define y as follows

y(1) is defined to be x

y(2) is defined to be xdot

Write out the derivative equations for y in terms of y. E.g.,

d( y(1) ) / dt = dx/dt = xdot = y(2)

d( y(2) ) / dt = d(xdot)/dt = xdotdot = ___________ <- this is what I asked you to solve for above

Just be sure to write that second equation in terms of y(1) and y(2). Once you have these equations figured out, you can write your code. Set initial values for y in your code and then write a loop to implement the Euler method.

Get started on this and then come back here to show us your progress and where you are stuck.

[Sander Z]

What about this code:

clear;

clc;

close('all');

% Mass-Spring System with Euler’s Method: (m((d^2)x)/(d(t^2)))+(c(dx/dt))+kx=0

% m = mass (kg)

% c = cushoning coefficient

% t = time (s)

% k = spring value (N/m)

% Model solution: ((mx'')+(cx')+(kx))=0

% m = 20;

% c = 40;

% k = 20;

% (20x'')+(10x')+(20x)=f

t = 0:0.1:15;

x0 = 1

x1 = 0

s12 = (-(c)+sqrt((c^2)-4*(m)*(k)))/(2*m); %find s21

s22 = (-(c)-sqrt((c^2)-4*(m)*(k)))/(2*m); %find s22

a12 = (-(s22)*x0+x1)/(s12-s22);

a22 = ((s12)*x0-x1)/(s12-s22);

xt = (a12)*exp(-t)+(t*(a22)*exp(-t))

plot(t,x);

title('For Critically Cushoning');

xlabel('X');

ylabel('y=f(x)');

grid;

What can I do to solve this? Am I completly wrong?

[John D'Errico]

Are you wrong? In the sense that you were apparently supposed to be solving it using Euler's method, yes you are wrong.

It looks like you are wanting to solve it using an analytical solution for x(t) though. Not sure if your solution is correct, but even if it is, you solve for x, as the variable xt. Then you plot x versus t, as opposed to xt versus t. You never actually defined any variable named x. So that plot must fail.

Anyway, IF you need to use Euler, then start at the beginning again, reading what James told you to do.

[Sander Z]

Forget about it, I though something else and went to the analytical mode for solving that, which is wrong.

I'll try as James told before.

Thanks for the tip.

[Sander Z]

I think the solution for the highest DE is like this one:

d(y(1))/dt = dx/dt = xdot = y(2)

d(y(2))/dt = d(xdot)/dt = xdotdot = (-1/m) * (c(dx/dt)+kx)

Right?

[James Tursa]

Almost. If you go with my suggestion of having a 2-element state vector y in your code, there won't be any dx/dt or x to use on the right hand side. That is where the definition of y(1) and y(2) come into play. Your second equation is correct, but you need to replace dx/dt and x by their equivalent y elements, in this case y(2) and y(1) respectively. So, taking it one more step:

d(y(2))/dt = d(xdot)/dt = xdotdot = (-1/m) * (c(dx/dt)+kx) = (-1/m)*(c*y(2) + k*y(1))

Thus, given the current state y, the two derivative equations then become:

dy(1) = y(2);

dy(2) = (-1/m)*(c*y(2) + k*y(1));

Here I am using a 2-element vector dy to contain the derivative of the 2-element current state vector y.

A very rough code outline would then be

Set initial constants m and k
Set initial y vector and initial time
Set dt value
Set n = number of iterations
for i=1:n
    dy(1) = (from above)
    dy(2) = (from above)
    (the next y) = (the current y) + dy * dt;
end

See if you can work on writing the code for the lines I have indicated above.

If you will be plotting the results, you will want to save all of those y values in a matrix. Let me know if you also need help with that (but first work on just getting the code above written).

[Sander Z]

Awsome! Now I get it what I should done, really good, ok, now I'll try another shoot in the dark, what about this code below with your function:

clear; 
clc; 
close('all');
% Mass-Spring System with Euler’s Method: (m((d^2)x)/(d(t^2)))+(c(dx/dt))+kx=0
Dt = 0.5; %response time [s]
m = 20; %mass [kg]
k = 20; %spring value [N/m]
c = 40; %cushoning value [Ns/m]
t0 = 0; %integration inicial time [s]
tf = 15; %integration final time [s]
t = t0:Dt:tf; %time vector
function dx=f(t,x)
    d(y1)/dt = dx/dt = xdot = y(2);
    d(y(2))/dt = d(xdot)/dt = xdotdot = (-1/m) * (c(dx/dt)+kx) = (-1/m)*(c*y(2) + k*y(1));
    dx = [y1; y2];
end
y0 = [1.0;0.0]; %inicial state y0 = [x;xp]
y = ode(y0, t0, t, f);

[James Tursa]

So, let's talk about your derivative function.

First, what you have written is not valid MATLAB syntax. What I had written above at first was just notional. Just use the dy(1) and dy(2) lines I had written.

Second, you have a function where you specify the input argument as x, but you use y in the body of the function. To keep it simple (and to match the notation that MATLAB uses in their doc), I would keep this as y and not use x. Also it is convenient to pass in the constants m, c, and k. So the function becomes

function dy = f(t,y,m,c,k)  % <-- changed x to y
    dy = [ y(2); (-1/m)*(c*y(2) + k*y(1)) ];
end

[John D'Errico]

Moved an answer to a comment (by Sander Z):

Got it, but now I'm a little confused, what would be the definition for the variable "x"?

The code must be like this below right?

clear; 
clc; 
close('all');
% Mass-Spring System with Euler’s Method: (m((d^2)x)/(d(t^2)))+(c(dx/dt))+kx=0
Dt = 0.5; %response time [s]
m = 20; %mass [kg]
k = 20; %spring value [N/m]
c = 40; %cushoning value [Ns/m]
t0 = 0; %integration inicial time [s]
tf = 15; %integration final time [s]
t = t0:Dt:tf; %time vector
y0 = [1.0;0.0]; %inicial state y0 = [x;xp]
n = 1; %number of iterations
y(1) = x;
y(2) = xdot;
for i=1:n
    dy(1) = y(2);
    dy(2) = (-1/m)*(c*y(2) + k*y(1));
 end

[James Tursa]

These lines:

n = 1; %number of iterations
y(1) = x;
y(2) = xdot;
for i=1:n
    dy(1) = y(2);
    dy(2) = (-1/m)*(c*y(2) + k*y(1));
 end

Should be something like this instead:

n = numel(t); %number of iterations
y = y0; % y0 contains the initial state
dy = zeros(2,1); % force dy to be a column vector
for i=1:n-1
    dy(1) = y(2);
    dy(2) = (-1/m)*(c*y(2) + k*y(1));
    y = y + dy * Dt; % you need to update y at each step using Euler method
 end

However, this will not store all the intermediate values of y ... it will simply overwrite y with the updated values. If you want to store the intermediate values (e.g., for plotting), you need to modify the above code to do so. Hint: One way to do it is to define y as a 2 x n matrix, where each column y(:,i) ( i.e., the elements y(1,i) and y(2,i) ) is the state at time t(i). See if you can implement that. You will need to change the syntax everywhere you use y to use two indexes.