System of 2nd order ODE with Euler.

Hi,

I had a system of 2 2nd order ODE.

I got to this point :

I need to find the approximate solutions of y2(t).

M1, M2, G, L1,L2 are variables given by the user.

These are the initials conditions which are given by the user also(i guess ?)

Im a bit lost in what should i do. I know how euler works but not with this type of system.

thanks !

[Jan]

The Euler integration worls exactly the same for systems of equations or single equations. Just write down the code and try it.

[Torsten]

Implicit or explicit Euler ?

[Guillaume Theret]

explicit ! Im writing some code, i'll send asap !

[Guillaume Theret]

after writing some code :

function main
  xinit = 0;
  xfinal = 3;
  h = .5;
  y0 = 1;
  [x,y] = euler_explicit(@fnc,xinit,xfinal,h,y0);
  plot(x,y(:,1));
end
function [x,y_e]=euler_explicit(f,xinit,xfinal,h,y0)
  x = xinit : h : xfinal;
  n = length(x);
  disp(n);
  y_e =zeros(1,n);
  disp(y_e)
  y_e(1) = y0;
  % compute y_e : euler explicit
  for i = 1:n - 1
    y_e(i+1) = y_e(i) + h *f(x(i),y_e(i));
  end
end
function dy = fnc(t,Y)
  L1 = 1;
  L2 = 2;
  M1 = 2;
  M2 = 3;   
  g = 1;
  K = 1/(L1*L2(M1+M2*sin(Y(1) - Y(2)).^2));
  Y4 = K*((M1+M2)*g*L1*sin(Y(1))*cos(Y(1)-Y(2)) - (M1+M2)*g*L1*sin(Y(2)) + (M1+M2)*L1^2*sin(Y(1) - Y(2))*Y(3)^2 + M2*L1*L2*sin(Y(1)-Y(2))*cos(Y(1)-Y(2))*Y(4)^2);
  Y3 = K*(-(M1+M2)*g*L2*sin(Y(1)) + M2*g*L2*sin(Y(2))*cos(Y(1)-Y(2)) - M2*L1*L2*sin(Y(1) - Y(2))*cos(Y(1)-Y(2))*Y(3)^2 - M2*L2^2*sin(Y(1)-Y(2))*Y(4)^2);
  dy = [Y(3),Y(4),Y3,Y4];
end

Im getting these erros. I check at my array and it looks like it has seven 0 which is normal.

I don't know if im going in the good direction or not to solve my equations :(

Thanks !

Answer 1

[Jan]

You got it almost. I've fixed a typo and expanded the Euler method to collect the output as matrix.

% function main
  xinit  = 0;
  xfinal = 3;
  h  = 0.05;
  y0 = [1, 0, 0, 0];  % As many elements as the system has
  [x, y] = euler_explicit(@fnc, xinit, xfinal, h, y0);
  plot(x, y);

% end

function [x, y] = euler_explicit(f, xinit, xfinal, h, y0)
  x = xinit : h : xfinal;
  n = length(x);
  y = zeros(n, numel(y0));
  y(1, :) = y0;
  for k = 1:n - 1
    y(k + 1, :) = y(k, :) + h * f(x(k), y(k, :));
  end
end

function dy = fnc(t,Y)
  L1 = 1;
  L2 = 2;
  M1 = 2;
  M2 = 3;   
  g = 1;
  K = 1 / (L1 * L2 * (M1 + M2*sin(Y(1) - Y(2)).^2));
  %                ^ was missing
  Y4 = K*((M1+M2)*g*L1*sin(Y(1))*cos(Y(1)-Y(2)) - (M1+M2)*g*L1*sin(Y(2)) + (M1+M2)*L1^2*sin(Y(1) - Y(2))*Y(3)^2 + M2*L1*L2*sin(Y(1)-Y(2))*cos(Y(1)-Y(2))*Y(4)^2);
  Y3 = K*(-(M1+M2)*g*L2*sin(Y(1)) + M2*g*L2*sin(Y(2))*cos(Y(1)-Y(2)) - M2*L1*L2*sin(Y(1) - Y(2))*cos(Y(1)-Y(2))*Y(3)^2 - M2*L2^2*sin(Y(1)-Y(2))*Y(4)^2);
  dy = [Y(3), Y(4), Y3, Y4];
end

[Guillaume Theret]

Thanks a lot !

Answer 2

[Torsten]

1. y0 must be a 4x1 vector, not a scalar.
2. ye = zeros(4,n) instead of ye=zeros(1,n)
3. ye(:,1) = y0 instead of ye(1) = y0
4. ye(:,i+1) = ye(:,i) + h*f(x(i),ye(:,i)) instead of the expression in your loop
5. dy = [Y(3);Y(4);Y3;Y4] instead of the row vector in your code

[Guillaume Theret]

Thanks. This is working. I had a syntax issue as mentionned by Jan