Eulers for two degree spring system

Hi there, I have to use basic Euler's to find and plot solutions to x1 and x2 over 0 < t < 30 where my step size is 0.001 with initial conditions: my conditions are x1(0) = 5, x2(0) = 1, v1(0) = 0,

Firstly I have multiplied the matrix attached out, subbed in v for x' and rearranged for x'' and renamed it as dv/dt where m = 1 c = 1 k = 2

 i have obtained: dv/dt = 1 - (4(x1) - 2(x2)) - (2(v1) - (v2))
                  dv/dt = 1 - (4(x2) - 2(x1)) - (2(v2) - (v1))

I then tried to implement Euler's for the first equation on matlab t=zeros(30000,0); y=zeros(30000,0);

t(0)=0; x(0)=5; v(0)=0;

for i=1:30000

    t(i+1)= t(i)+0.001;
    x(i+1)= x(i)+0.001.*(1-(4*x(i)-2*x(i)) -(2*v(i)-v(i));

end

plot(t,x)

I'm getting an error in line 11, unbalances parenthesis and I'm not sure I'm really implementing the scheme correctly so any help would be greatly appreciated!

Answer 1

[James Tursa]

You don't have enough equations in your code. This is a 4th order system (Two variables each as a 2nd order ODE gives 2x2 = 4 orders), so you will need 4 1st order equations. Using the nomenclature in the MATLAB doc for the ODE methods, let's pick y to be the state vector of your system (this will make it easier to compare your code with the ODE examples in the doc). Since you have a 4th order system, y will have 4 elements. Let's pick the elements to be the following:

y(1) = x1
y(2) = x2
y(3) = d(x1)/dt = x1dot = v1
y(4) = d(x2)/dt = x2dot = v2

The derivative of your state vector would then be

d(y(1))/dt = d(x1)/dt = y(3)
d(y(2))/dt = d(x2)/dt = y(4)
d(y(3))/dt = something
d(y(4))/dt = something

For those last two derivative equations, the "something" on the right hand side will be what you get by solving your two posted 2nd order equations for x1dotdot and x2dotdot and then substituting the appropriate y elements for the x1, x2, x1dot, and x2dot variables.

Then for your Euler iterations, at each step you will be propagating this 4 element vector. To do that, one way is to redefine our y vector as a matrix where each column is a y vector at a given time point. E.g.,

y(1,i) = x1 at time t(i)
y(2,i) = x2 at time t(i)
y(3,i) = x1dot at time t(i)
y(4,i) = x2dot at time t(i)

Then using this definition and the equations discussed above your MATLAB code would look something like this:

t = zeros(30001); % pre-allocate the time vector
y = zeros(4,30001); % pre-allocate the state vector
h = 0.001; % The stepsize
t(1) = 0; % Initial time (needs to start at index 1, not index 0)
y(1,1) = 5; % Initial x1 at time t(1) is contained in first column at y(1,1) 
y(2,1) = 1; % Initial x2 at time t(1) is contained in first column at y(2,1) 
y(3,1) = 0; % Initial v1 at time t(1) is contained in first column at y(3,1) 
y(4,1) = 0; % Initial v2 at time t(1) is contained in first column at y(4,1) 
for i=1:30000 % Note that the i=1 spot is already filled with initial conditions
    t(i+1)= t(i) + h;
    y(1,i+1) = y(1,i) + h * (     );
    y(2,i+1) = y(2,i) + h * (     ); 
    y(3,i+1) = y(3,i) + h * (     );
    y(4,i+1) = y(4,i) + h * (     );
end
plot(t,y(1,:),t,y(2,:)); % y(1,:) is the x1 solution, and y(2,:) is the x2 solution

What you need to do is fill in the blanks inside the parentheses in the code above. Everything inside those parentheses will be an expression involving y(1,i), y(2,i), y(3,i), and y(4,i). The first two should be easy based on the above discussion. The last two will involve your solution to solving your two posted 2nd order equations. Give it a try, and if you get stuck go ahead and post how far you got and what you are stuck on.

[gary stewart]

thanks for your answer but I'm still struggling within the brackets;

y(1,i+1) = y(1,i) + h * ((4*y(1,l)-2*y(2,l))-(2*y(3,l)-(4,l));
y(2,i+1) = y(2,i) + h * ((4*y(2,l)-2*y(1,l))-(2*y(4,l)-(3,l));

is this correct? if so i have no idea how I would then use this to find equations for y3 and y4

thanks again!

[gary stewart]

      y(1,i+1) = y(1,i) + h * (y(1,l)); as this equals x1?
      y(2,i+1) = y(2,i) + h * (y(2,l)); as this equals x2?
      y(3,i+1) = y(3,i) + h * ((4*y(1,l)-2*y(2,l))-(2*y(3,l)-(4,l)); this equals v1?
      y(4,i+1) = y(4,i) + h * ((4*y(2,l)-2*y(1,l))-(2*y(4,l)-(3,l)); this equals v2?

[James Tursa]

Tackle the first two equations first. From our my definition of y above and adding in the column index i we have this:

d(y(1,i))/dt = d(x1_i)/dt = y(3,i)
d(y(2,i))/dt = d(x2_i)/dt = y(4,i)

So that tells you what to put in the parentheses for your first two equations. You get this:

      y(1,i+1) = y(1,i) + h * (y(3,i));
      y(2,i+1) = y(2,i) + h * (y(4,i));

For the last two equations, you use the solution for the two 2nd order equations you did above, but again write it in terms of y(1,i), y(2,i), y(3,i) and y(4,i). E.g.,

d(y(3,i))/dt = d(x1dot at time t(i))/dt = x1dotdot at time t(i)
d(y(4,i))/dt = d(x2dot at time t(i))/dt = x2dotdot at time t(i)

Then the expressions for x1dotdot (which is v1dot) and x2dotdot (which is v2dot) is from your solution to the 2nd order equations. You just need to plug in y(1,i), y(2,i), y(3,i) and y(4,i) for the appropriate x1, x2, v1, v2 elements.