code to perform euler method not working properly

Question

0 个投票

Hi all

I'm writing a code to solve a system of 8 first order differential equations(which I transformed from 4 second order equations) using the Euler method. The solution should show 4 lines, the first three show the position of the front tire, rear tire, and chassis of a motorbike after an initial bump is hit by the bike. the fourth line gives the angle of the chassis in respect to the horizontal. The code solves the system and for a small time it looks like it could possibly be the correct solution but after a few seconds its clear that the solution does not make sense... If someone could take a look at the code I wrote and confirm that it SHOULD work that would be greatly appreciated, if that is the case then I must have made a mistake deriving the equations at some point.

Thanks in advance.

function dx = projectmatrix2(t,u1)
k1 = 10;
k2 = 10;
k3 = 10;
m1 = 2.5;
m2 = 2.5;
m3 = 3;
c1 = 15;
c2 = 15;
r = 1;
I = 22;
L = 4;
dx = zeros(8,1);
dx(1) = u1(2);
dx(2) = -(k1*u1(1)-c1*(u1(2)-u1(6))+k1*(L+u1(1)-u1(5)))/m1;
dx(3) = u1(4);
dx(4) = -(k2*u1(5)-c2*(u1(4)-u1(6))+k2*(L+u1(3)-u1(5)))/m2;
dx(5) = u1(6);
dx(6) = (c1*(u1(2)-u1(6))+c2*(u1(4)-u1(6))-k1*(L+u1(1)-u1(5))-k2*(L+u1(3)-u1(5)))/m3;
dx(7) = u1(8);
dx(8) = (r*(c1*(u1(2)-u1(6))+c2*(u1(4)-u1(6))-k1*(L+u1(1)-u1(5))-k2*(L+u1(3)-u1(5))))/I;
end

The m. file above is saved and then referenced in the file below.

%Euler method
Tsim = 10;
h = .1;
n = Tsim/h;
x = [1;0;0;0;3;0;0;0];
t = zeros(1,n);
for i=1:n-1
   t(i+1) = t(i) + h;
   x(:,i+1) = x(:,i) + h * projectmatrix2(t(i),x(:,i));
end
plot(t,x(1,:),'b*-')
hold on
plot(t,x(3,:),'g*-')
plot(t,x(5,:),'r*-')   
plot(t,x(7,:),'k*-')
hold off