I am currently working on some system of ODEs enter image description here.

These four odes are equation of motion for a physical system, photo linked below:

Scroll down to the bottom of this page, the physical system in motion and interactive java applet: https://www.myphysicslab.com/pendulum/cart-pendulum-en.html.

x represents displacement of the main block q represents the angle theta of the pendulum from the vertical position v represents velocity of the main block w represents the angular velocity omega of the pendulum.

I coded myself a matlab file to draw a graph of x vs t(time), however, it did not work and produced extremely large outputs. I wonder if anyone could help me with the code (since I am a newbie) and give me some advices.

Btw, I chose to use this numerical method because the equations give implicit solutions when solved normally (nonlinear).

% x'= v

% q'= w

% v' = f(x,q,v,w) = (pR(w^2)sin(q)+pg*sin(q)cos(q)-kx-c*v+(b/R)wcos(q))/((m+p*(sin(q))^2))

% w' = g(x,q,v,w) = (-pR(w^2)*sin(q)*cos(q)-(p+m)gsin(q)+kxcos(q)+cvcos(q)-(1+M/q)*(b/R)w)/(R(M+p*(sin(q))^2))

p = 0.05; % mass of pendulum mass

m = 0.198; % mass of large mass

R = 0.1; % length of pendulum rod

g = 9.80665; % gravitational acceleration

k = 15; % large mass spring constant

c = 0.124; % large mass damping coefficient

b = 0.00543 % damping torque coefficient for pendulum mass

% Equations

fX = @(t,V) V;

fQ = @(t,W) W;

fV = @(t,X,Q,V,W) (pR(W^2)sin(Q)+pg*sin(Q)cos(Q)-kX-c*V+(b/R)Wcos(Q))/((m+p*(sin(Q))^2));

fW = @(t,X,Q,V,W) (-pR(W^2)*sin(Q)*cos(Q)-(p+m)gsin(Q)+kXcos(Q)+cVcos(Q)-(1+m/p)*(b/R)W)/(R(m+p*(sin(Q))^2));

% Initial conditions

t(1) = 0;

X(1) = 0.1;

Q(1) = 0;

V(1) = 0;

W(1) = 0;

% Step size and number of steps h = 0.0001; time_total = 20; n = ceil(time_total/h);

% RK4 loop for i = 1:n

t(i+1) = t(i) + h; % time increment

k1X = fX(t(i),V(i));
k1Q = fQ(t(i),W(i));
k1V = fV(t(i),X(i),Q(i),V(i),W(i));
k1W = fW(t(i),X(i),Q(i),V(i),W(i));

k2X = fX(t(i)+h/2,V(i)+h/2+k1V);
k2Q = fQ(t(i)+h/2,W(i)+h/2+k1W);
k2V = fV(t(i)+h/2,X(i)+h/2+k1X,Q(i)+h/2+k1Q,V(i)+h/2+k1V,W(i)+h/2+k1W);
k2W = fW(t(i)+h/2,X(i)+h/2+k1X,Q(i)+h/2+k1Q,V(i)+h/2+k1V,W(i)+h/2+k1W);

k3X = fX(t(i)+h/2,V(i)+h/2+k2V);
k3Q = fQ(t(i)+h/2,W(i)+h/2+k2W);
k3V = fV(t(i)+h/2,X(i)+h/2+k2X,Q(i)+h/2+k2Q,V(i)+h/2+k2V,W(i)+h/2+k2W);
k3W = fW(t(i)+h/2,X(i)+h/2+k2X,Q(i)+h/2+k2Q,V(i)+h/2+k2V,W(i)+h/2+k2W);

k4X = fX(t(i)+h,V(i)+h+k3V);
k4Q = fQ(t(i)+h,W(i)+h+k3W);
k4V = fV(t(i)+h,X(i)+h+k3X,Q(i)+h+k3Q,V(i)+h+k3V,W(i)+h+k3W);
k4W = fW(t(i)+h,X(i)+h+k3X,Q(i)+h+k3Q,V(i)+h+k3V,W(i)+h+k3W);

X(i+1) = X(i)+(h/6)*(k1X+2*k2X+2*k3X+k4X);
Q(i+1) = Q(i)+(h/6)*(k1Q+2*k2Q+2*k3Q+k4Q);
V(i+1) = V(i)+(h/6)*(k1V+2*k2V+2*k3V+k4V);
W(i+1) = W(i)+(h/6)*(k1W+2*k2W+2*k3W+k4W);



  • $\begingroup$ Please indent all of your code so that it formats as code block, read the guide on the markup (markdown?). $\endgroup$ – LutzL Aug 10 '17 at 9:43

One obvious problem is that

k2X = fX(t(i)+h/2,V(i)+h/2+k1V);

should read

k2X = fX(t(i)+h/2,V(i)+h/2*k1V);

multiplication, not addition of $h/2$.

The general critique is that you are using the Matrix Laboratory, so in keeping with its philosophy, use vector valued functions and the built-in vector operations. Also use the provided ODE integrators and make your interfaces similar to them so that switching back and forth is easy.

  • $\begingroup$ Thanks, detail does matter! $\endgroup$ – Harry Yang Aug 11 '17 at 0:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.