Runge-Kutta Method solving four ODES

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

%RK1
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));

%RK1
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);

%RK1
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);

%RK1
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);

%Summing
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);


end

plot(t,X);

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

One obvious problem is that

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


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.

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