code to perform euler method not working properly

조회 수: 9 (최근 30일)
Hal Reigi
Hal Reigi 2017년 4월 18일
답변: Jan 2017년 4월 18일
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

채택된 답변

Jan
Jan 2017년 4월 18일
The shown behaviour occurres for h=0.01 and h=0.001 also. I assume, the formula in projectedmatrix2 is wrong.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by