numerical methodsordinary differential equationsrunge-kutta-methods
I have just started getting into ODEs, and have come across the Runge-Kutta method for numerically solving them. However, in playing around with them to model hypothetical situations, I came across the equations:
$$\begin{aligned} \dot x &= x – x^2 – y\\ \dot y &= xy – y^2 \end{aligned}$$
I was trying to think of how one would use Runge Kutta methods to do this, and I couldn't figure it out. What is this type of differential equation called, and how can I simulate/solve it?
Since you have not given the set of differential equations, this answer cannot respond exactly. However, as indicated in the comments, the system of equations is implicit. These can be written as $F(t,y,y')=0$, where $$y(t)=[y_1(t),y_2(t),\ldots,y_n(t)]^T \text{ and } y'(t)=[y'_1(t),y'_2(t),\ldots,y'_n(t)]^T.$$This system can be solved numerically in MATLAB using ode$15$i.
As an example, suppose we want to solve the implicit system of ODEs $$ F(t,y,y')=\begin{bmatrix}\cos(y'_1)+y'_2-y_1\\t\sin(y'_2)+y'_1-y_2^2\end{bmatrix} $$ with $y(0)=[0,0]^T$ for $0\leq t\leq T$.
ode$15$i requires consistent initial conditions for both $y$ and $y'$. Given that we do not have initial conditions, they can be approximated using the decic command.
This can be done in MATLAB as follows.
T = 5;
y0 = [0;0];
fun = @(t,y,yp) [cos(yp(1))+yp(2)-y(1);t*sin(yp(2))+yp(1)-y(2)^2];
[y0,yp0] = decic(fun,0,y0,[1;1],[0;0],[0;0]);%get consistent initial conditions
[tsol,ysol] = ode15i(fun,[0,T],y0,yp0); % solve the system
plot(tsol,ysol(:,1),'b',tsol,ysol(:,2),'r') % plot the solution
Edit: Solution to include equations
System
The system of DEs to be solved is as follows: $$ \begin{align} \theta '' &= \frac{\theta' \cos \theta}{R} + \frac{R'\sin \theta}{R} + \frac{\nu\sin\theta}{\kappa R} + \frac{\eta\cos\theta}{\kappa R} \tag{1}\\ \nu' &= \frac{\kappa}{2} \left( \theta' + \frac{\sin\theta}{R}-c_0 \right)\left( \theta'-\frac{\sin\theta}{R}-c_0\right) + \lambda\tag{2}\\ \eta' & =0\tag{3}\\ R' &= \cos \theta \tag{4}\\ z' &=- \sin \theta \tag{5}\end{align} $$ No initial conditions are given, though it is stated that all zero initial conditions may be appropriate.
Solution
Noting $(3)$, the solution of $\eta'=0$ must be $\eta(t)=c$ for some constant $c$. Since $\eta(0)=0$, we must have $c=0$, giving $\eta(t)=0$. As such, $(3)$ is not needed, since $\eta(t)$ can be obtained analytically. I will assume, however, throughout the remainder of the solution below, that $\eta(t)$ need not be a zero constant, but simply any constant, as given by the initial conditions.
Furthermore, $(1)$ can be reduced to first order by introducing the auxiliary variable $A=\theta'$. Then $(1)$ becomes $$ A'=\frac{A \cos \theta}{R} + \frac{R'\sin \theta}{R} + \frac{\nu\sin\theta}{\kappa R} + \frac{\eta\cos\theta}{\kappa R} \tag{6} $$ with the additional equation $$ \theta'=A\tag{7}. $$ The system then required to be solved is $(2),\,(4),\,(5),\,(6)\text{ and }(7)$.
In MATLAB, an m-file for the ode function should be created similar to as follows
function F = odefun(t,y,yp,parms)
% evaluate the ode functions with the given parameters
% y and yp are ordered as [theta,nu,R,z,A]
% parms is a structure, containing eta, kappa, lambda, c0
% extract components from parms structure
eta = parms.eta;
kappa = parms.kappa;
lambda = parms.lambda;
c0 = parms.c0;
% extract components from y and yp for more clear use
theta = y(1); thetap = yp(1);
nu = y(2); nup = yp(2);
R = y(3); Rp = yp(3);
z = y(4); zp = yp(4);
A = y(5); Ap = yp(5);
% allocate memory for return variable
F = zeros(5,1);
% precompute some quantities
sint = sin(theta);
cost = cos(theta);
% evaluate equation (2)
F(1) = nup - (kappa/2*(thetap+sint/R-c0)*(thetap-sint/R-c0)+lambda); %****
% evaluate equation (4)
F(2) = Rp - cost;
% evaluate equation (5)
F(3) = zp + sint;
% evaluate equation (6)
F(4) = Ap - (A*cost + Rp*sint + (nu*sint+eta*cost)/kappa)/R; %****
% evaluate equation (7)
F(5) = thetap - A;
Note: Both equations marked by %**** $\left(\text{i.e. $(2)$ and $(6)$}\right)$ will have a divide by zero if $R(0)=0$.
From my previous smaller example, you should be able to make this function work with the ode15i(...) and decic(...) commands.
There are many Runge–Kutta methods. The one you have described is (probably) the most popular and widely used one. I am not going to show you how to derive this particular method – instead I will derive the general formula for the explicit second-order Runge–Kutta methods and you can generalise the ideas.
In what follows we will need the following Taylor series expansions.
\begin{align} \tag1\label{eq1} & x(t+\Delta t) = x(t) + (\Delta t)x'(t) + \frac{(\Delta t)^{2}}{2!}x''(t) + \text{higher order terms}. \\ \tag2\label{eq2} & f(t + \Delta t, x + \Delta x) = f(t,x) + (\Delta t)f_{t}(t,x) + (\Delta x)f_{x}(t,x) + \text{higher order terms}. \end{align}
We are interested in the following ODE:
$$x'(t) = f(t,x(t)).$$
The value $x(t)$ is known and $x(t+h)$ is desired. We can solve this ODE by integrating: $$x(t+h) = x(t) + \int_{t}^{t+h}f(\tau,x(\tau))\text{d}\tau.$$ Unfortunately, actually doing that integration exactly is often quite hard (or impossible), so we approximate it using quadrature: $$x(t+h) \approx x(t) + h\sum_{i=1}^{N}\omega_{i}f(t + \nu_{i}h,x(t + \nu_{i}h)).$$ The accuracy of the quadrature depends on the number of terms in the summation (the order of the Runge–Kutta method), the weights, $\omega_{i}$, and the position of the nodes, $\nu_{i}$.
Even this quadrature can be quite difficult to calculate, since on the right-hand side we need $x(t + \nu_{i}h)$, which we don’t know. We get around this problem in the following manner:
Let $\nu_1=0$, which makes the first term in the quadrature $K_{1} := hf(t,x(t)).$ This we do know and we can also use it to approximate $x(t + \nu_{2}h)$ by writing the second term in the quadrature as $K_{2} := hf(t+\alpha h,x(t) + \beta K_{1})$.
With this, the quadrature formula is: $$\tag3\label{eq3} x(t+h) = x(t) + \omega_{1}K_{1} + \omega_{2}K_{2}.$$
Some notes:
In equation \eqref{eq3}, we substitute in the Taylor series expansion \eqref{eq1} on the left-hand side:
$$x(t) + h x'(t) + \frac{h^{2}}{2!}x''(t) + \mathcal{O}\left(h^{3}\right) = x(t) + \omega_{1}K_{1} + \omega_{2}K_{2}.$$
Since $x' = f$ and thus $x'' = f_{t} + ff_{x}$ (suppressing arguments for ease of notation), we have:
$$hf + \frac{h^{2}}{2}(f_{t} + ff_{x}) + \mathcal{O}\left(h^{3}\right) = \omega_{1}K_{1} + \omega_{2}K_{2}.$$
Now substitute for $K_{1}$ and $K_{2}$:
$$hf + \frac{h^{2}}{2}(f_{t} + ff_{x}) + \mathcal{O}\left(h^{3}\right) = \omega_{1}hf + \omega_{2}hf(t + \alpha h, x + \beta K_{1}).$$
Taylor-expand on the right-hand side using \eqref{eq2}:
$$hf + \frac{h^{2}}{2}(f_{t} + ff_{x}) + \mathcal{O}\left(h^{3}\right) = \omega_{1}hf + \omega_{2}(hf +\alpha h^{2}f_{t} + \beta h^{2} ff_{x}) + \mathcal{O}\left(h^{3}\right).$$
Thus the Runge–Kutta method will agree with the Taylor series approximation to $\mathcal{O}\left(h^{3}\right)$ if we choose:
$$\omega_{1} + \omega_{2} = 1,$$ $$\alpha \omega_{2} = \frac{1}{2},$$ $$\beta \omega_{2} = \frac{1}{2}.$$
The canonical choice for the second-order Runge–Kutta methods is $\alpha = \beta = 1$ and $\omega_{1} = \omega_{2} = 1/2.$
The same procedure can be used to find constraints on the parameters of the fourth-order Runge–Kutta methods. The canonical choice in that case is the method you described in your question.
Best Answer
The Runge-Kutta formulas for a system of differential equations are really the same as for a single equation, it's just that your dependent variable is a vector rather than a scalar. Write your system as $$\dfrac{dX}{dt} = F(t, X(t))$$ where $X = (x, y)$. If you're using the classical fourth-order Runge-Kutta with step size $h$, your iteration is $$ \eqalign{K_1 &= h F(t_n, X_n)\cr K_2 &= h F(t_n + h/2, X_n + K_1 /2)\cr K_3 &= h F(t_n + h/2, X_n + K_2 /2)\cr K_4 &= h F(t_n + h, X_n + K_3)\cr t_{n+1} &= t_n + h\cr X_{n+1} &= X_n + (K_1 + 2 K_2 + 2 K_3 + K_4)/6\cr}$$