MATLAB Answers

System of 2nd order ODE with Euler.

13 views (last 30 days)
Guillaume Theret
Guillaume Theret on 27 May 2021
Commented: Guillaume Theret on 27 May 2021
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 !
  4 Comments
Guillaume Theret
Guillaume Theret on 27 May 2021
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 !

Sign in to comment.

Accepted Answer

Jan
Jan on 27 May 2021
Edited: Jan on 27 May 2021
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

More Answers (1)

Torsten
Torsten on 27 May 2021
  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
  1 Comment
Guillaume Theret
Guillaume Theret on 27 May 2021
Thanks. This is working. I had a syntax issue as mentionned by Jan

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!