Numerical solution of an ODE system of equations using RK4

Question

I have given an assignment to find the solution to the ODE system of equations as follow:

$$\begin{cases} x_1' = x_1 + x_2 \\ x_2' = -3x_1 -10x_2 + x_2 ^2\end{cases}$$

With initial conditions: $x_1(0) = 1, \quad x_2(0)=10$ on interval $[0,10]$

The above system of questions differs in two aspects from the normal questions of my textbook so I am confused and don't know how to program the algorithm in Matlab:

1. It has two intermingled equations. In contrary to textbook samples which has only one in form of: $x'=f(t,x)$
2. $x_1'$ and $x_2'$ don't have any $t$ variable! which makes me even more confused.

I appreciate it if you elaborate on the above two points. Thanks.

------------------------------------------Part2-----------------------------------------

Ok with above clarifications I started to program in Matlab, after running my code, $x_1(t)$ and $x_2(t)$ seems to change (probably becomes unstable) after changing the step number which changes h respectively.

Here is the code:

function [sol_x1, sol_x2] = rk4 (m)
    a = 0;      % interval = [a, b] = [0, 10]
    b = 10;
    h = (b-a) / m; % m is  number of steps

    Sx1 = zeros (1,m+1);      % Solutions will be saved here
    Sx2 = zeros (1,m+1);

    Sx1(1) = 1;    % initial values
    Sx2(1) = 10;

    for j = 1:m
        x1 = Sx1(j);
        x2 = Sx2(j);

        fprintf('x1 = %i, x2 = %i \n', x1, x2);

        x1_k1 = h * f1(x1, x2);
        x2_k1 = h * f2(x1, x2);

        x1_k2 = h * f1(x1 + x1_k1/2, x2 + x2_k1/2);
        x2_k2 = h * f2(x1 + x1_k1/2, x2 + x2_k1/2);

        x1_k3 = h * f1(x1 + x1_k2/2, x2 + x2_k2/2);
        x2_k3 = h * f2(x1 + x1_k2/2, x2 + x2_k2/2);

        x1_k4 = h * f1(x1 + x1_k3, x2 + x2_k3);
        x2_k4 = h * f2(x1 + x1_k3, x2 + x2_k3);


        next_x1 = x1 + (x1_k1 + (2 * x1_k2) + (2 * x1_k3) + x1_k4) / 6;
        next_x2 = x2 + (x2_k1 + (2 * x2_k2) + (2 * x2_k3) + x2_k4) / 6;

        % no need for t = t + h because x1 and x2 do not depend on t

        %Save the solutions
        fprintf('%i\n', h);
        %display(h)
        Sx1(j+1) = next_x1;
        Sx2(j+1) = next_x2;
    end    

    sol_x1 = Sx1;
    sol_x2 = Sx2;
end

function [y] = f1 (x1, x2) 
    y = x1 + x2;
end

function [y] = f2 (x1, x2) 
    y = (-3 * x1) + (-10 * x2) + (x2 ^ 2);
end

Am I on right path?

$x1(t)$ and $x2(t)$ shoots up to very large numbers pretty quickly. I think something is wrong with my algorithm implementation.

Answer 1

Additionally, you have a problem with the line

h = h + step

while you use h as step size. Since originally step = h you get a rapidly growing step size where it should be constant. What you probably meant was

t = t + h

there is no need for the extra step variable. Also, consider using h=(b-a)/m for greater flexibility, using a single point of constant definition.

---

Everything else is looking good. You can always test for reliability of the code by varying the step size. The numerical results should stay stable with small variations.

Answer 2

Note that the notation $f(x,t)$ is for completeness, in the sense that it is also meant to cover cases where there is no explicit $t$ dependence.

As for the other part of your question, you can write

\begin{align} x_1' & = f_1(x_1,x_2,t) = \cdots \\ x_2' & = f_2(x_1,x_2,t) = \cdots \end{align}

If you want me to see if I can find some examples online, or write some pseudo-code, I can do that for you.

Answer 3

You can think of this system of differential equations as ${\bf x}' = f({\bf x},t)$ where ${\bf x}$ is the vector $\pmatrix{x_1\cr x_2\cr}$, and $$ f\left(\pmatrix{x_1\cr x_2\cr},t\right) = \pmatrix{x_1 + x_2\cr −3 x_1−10 x_2+x_2^2}$$ The fact that this doesn't happen to depend on $t$ makes no difference. The fact that $\bf x$ is a vector rather than a scalar also doesn't matter: the RK4 formulas are exactly the same as in the scalar case.