Solving ODEs in MATLAB using the Runga - Kutta Method

Question

I am required to solve these particular ODEs using numerical methods in MATLAB. The ODEs essentially model the fall of a body of mass m, connected to a piece of elastic with spring constant k. The solutions to these ODEs are to represent the body's position and velocity at discrete positions in time.

The parameters for the ODEs are,

 H = 74             
 D = 31              
 c = 0.9          
 m = 80             
 L = 25             
 k = 90              
 g = 9.8            
 C = c/m             
 K = k/m            
 T = 60             
 n = 10000           

I've implemented the following two methods; Euler and Fourth Order Runge - Kutta to approximate the solutions over the interval [0, 60].

---

Here is my Runge Kutta function,

function [t, y, v, h] = rk4_approx(T, n, g, C, K, L)

    %% calculates interval width h for each iteration
    h = (T / n);

    %% creates time array t
    t = 0:h:T;

    %% initialises arrays y and v to hold solutions
    y = zeros(1,n+1);
    v = zeros(1,n+1);

    %% functions
    z = @(v) v;
    q = @(v, y) (g - C*abs(v)*v - max(0, K*(y - L)));

    %% initial values
    v(1) = 0;
    y(1) = 0;

    %% performs iterations
    for j = 1:n

        %% jumper's position at each time-step
        r1 = h*z(v(j));
        r2 = h*z(v(j) + 0.5*h);
        r3 = h*z(v(j) + 0.5*h);
        r4 = h*z(v(j) + h);
        y(j+1) = y(j) + (1/6)*(r1 + 2*r2 + 2*r3 + r4); %position solution

        %% jumper's velocity at each time-step
        k1 = h*q(v(j), y(j));
        k2 = h*q(v(j) + 0.5*h, y(j) + 0.5*k1);
        k3 = h*q(v(j) + 0.5*h, y(j) + 0.5*k2);
        k4 = h*q(v(j) + h, y(j) + k3);
        v(j+1) = v(j) + (1/6)*(k1 + 2*k2 + 2*k3 + k4); %velocity solution

    end

end

---

Here is my Euler function,

function [t, y, v, h] = euler_approx(T, n, g, C, K, L)

    % calculates interval width h
    h = T / n;

    % creates time array t
    t = 0:h:T;

    % initialise solution arrays y and v
    y = zeros(1,n+1);
    v = zeros(1,n+1);

    % perform iterations
    for j = 1:n
        y(j+1) = y(j) + h*v(j);
        v(j+1) = v(j) + h*(g - C*abs(v(j))*v(j) - max(0, K*(y(j) - L)));
    end

end

---

However, after varying the parameter 'n' (where n is the number of 'steps' in the iteration) it appears the Euler solution for the position of the body converges to the maximum value of approximately y = 50 faster then the Runge - Kutta solution does. Since this ODE does not have a closed solution I have nothing to compare my answer to. I suspect the answer to be y = 50 though.

Therefore, I'm doubting my answer.

Is my code for the Runge - Kutta solution incorrect? Should it not converge faster than the Euler solution?

Sorry for my poor formatting.

Answer 1

The Runge-Kutta integration is incorrect.

It is vital to appreciate the difference between independent and dependent (also called state and a host of other names) variables. Time is the independent variable in this problem. Given a time, you can provide a height and a velocity; the reverse is not uniquely true. When you read a Runge-Kutta formula, such as the one provided by Wikipedia, t is the independent variable and y is vector of dependent variables. Also, when performing time integration of systems of equations (here we have a system of two equations), it is very important to keep track of which right-hand side belongs to which equation if you are going to perform the march element-wise, which I will do for simplicity.

All that said, the problem with the current RK integrator is two-fold

• v is being stepped as if it were t; this is incorrect. Both v and y are stepped similarly.
• y should be stepped with the r variables since the r variables come from y's right-hand side equation z. Similarly, v is stepped with the k variables.

The updated core of the integrator is thus:

r1 = h*z(v(j));
k1 = h*q(v(j), y(j));

r2 = h*z(v(j) + 0.5*k1);
k2 = h*q(v(j) + 0.5*k1, y(j) + 0.5*r1);

r3 = h*z(v(j) + 0.5*k2);
k3 = h*q(v(j) + 0.5*k2, y(j) + 0.5*r2);

r4 = h*z(v(j) + k3);
k4 = h*q(v(j) + k3, y(j) + r3);

y(j+1) = y(j) + (1/6)*(r1 + 2*r2 + 2*r3 + r4); %position solution
v(j+1) = v(j) + (1/6)*(k1 + 2*k2 + 2*k3 + k4); %velocity solution

Notice how both v and y are updated in a similar fashion and, therefore, are required to be updated in lock-step with one another. This form of the integrator will give far better performance than Euler.

Finally, if in doubt in the future about a solution you don't know, always remember you have the MATLAB ODE suite at your disposal, and a quick call to the extensively vetted and very robust ode45 can relieve a lot of concerns. I actually used this call

[t45,w45] = ode45(@(t,w) [z(w(2));q(w(2),w(1))],linspace(0,T,200).',[0;0]);

to check my work.