Bring ODE into a suitable form to solve it with Runge-Kutta steps

303 Views Asked by At

Can anyone please help me understand, how I should bring this ODE

y'' + y = sin(t)

with initial conditions y(0) = 100, y'(0) = 5

into a Runge-Kutta-Form?

I tried to solve this equation, the solution of the homogeneous equation is $$y_{\rm hom} = a + b\,e^{-x}$$ for some constants $a,b$. The particular solution is $y = 0$ (?)

Thank you very much

2

There are 2 best solutions below

0
On BEST ANSWER

In order to apply Runge-Kutta methods you need to write the given initial-value problem in the form \begin{eqnarray} \boldsymbol{\dot{y}} &=& \boldsymbol{f}(t,\boldsymbol{y}),\\ \boldsymbol{y}(t_0) &=& \boldsymbol{y}_0, \end{eqnarray} for some time-dependent, vector-valued function $\boldsymbol{y}(t)$. The function $\boldsymbol{f}$, the initial time $t_0$ and the initial value $\boldsymbol{y}_0$ need to be known. Then you can apply some s-stage Runge-Kutta method with parameters $a_{ij}, b_i, c_i$ and with step size $h > 0$: \begin{eqnarray} \boldsymbol{k}_i &=& \boldsymbol{f}\left(t_{n-1} + c_i h, \boldsymbol{y}_{n-1} + h \sum_{j=1}^s a_{ij} \boldsymbol{k}_j \right), \quad i=1,2,\dots,s,\\ \boldsymbol{y}_n &=& \boldsymbol{y}_{n-1} + h \sum_{i=1}^s b_i \boldsymbol{k}_i, \end{eqnarray} for $n=1,2,\dots$, to obtain approximations $\boldsymbol{y}_n \simeq \boldsymbol{y}(t_n)$, where $t_n = t_0 + n h$.

In order to obtain the required form of your problem you define the vector-valued function $\boldsymbol{y} := (y,\dot{y})^{\top}$, and you write using the given ODE: \begin{equation} \boldsymbol{\dot{y}} = \left( \begin{array}{c} \dot{y}\\ \ddot{y} \end{array} \right) = \left( \begin{array}{c} \dot{y}\\ - y + \sin(t) \end{array} \right) =: \boldsymbol{f}(t,\boldsymbol{y}). \end{equation} The initial point $(t_0, \boldsymbol{y}_0)$ is obtained from the given initial conditions: with $t_0 := 0$ you write \begin{equation} \boldsymbol{y}(0) = \left( \begin{array}{c} y(0)\\ \dot{y}(0) \end{array} \right) = \left( \begin{array}{c} 100\\ 5 \end{array} \right) =: \boldsymbol{y}_0. \end{equation} Now you have all the ingredients required in order to solve the problem numerically using a Runge-Kutta method.

This derivation also explains the code given by Max Herrmann, and in his explanation he used some more complicated version of an (embedded) Runge-Kutta method which features automatic step size selection.

0
On

For solving with Octave, enter the following into the command prompt:

f = @(t,y) [y(2); -y(1) + sin(t)];
[t,y] = ode45 (f, [0, 2*pi], [100, 5]);
plot(t,y)

The result should look something like this:

enter image description here

The blue line represents $y(t)$, the red one its time derivative.


The numeric solver uses a Runge-Kutta-based procedure (Dormand-Prince) to evaluate it numerically, by solving

$$ \begin{eqnarray} y_1^{(n+1)} & = & y_1^{(n)} + h\sum_{j=1}^{s}b_{j}k_{1j}^{(n)} \\ y_2^{(n+1)} & = & y_2^{(n)} + h\sum_{j=1}^{s}b_{j}k_{2j}^{(n)}, \end{eqnarray} $$

where $h$ is the time increment, $s$ and $b$ are integration-method-specifics, and $k_{ij}$ is evaluated with the right-hand side as detailed below.


The classical Runge-Kutta method, aka RK4, uses the parameters $s=4$, $b_1 = b_4 = \frac{1}{6}$, $b_2 = b_3 = \frac{1}{3}$, and

$$ \begin{eqnarray} k_{11}^{(n)} & = & y_2^{(n)} \\ k_{21}^{(n)} & = & -y_1^{(n)} + \sin(t_n) \\ k_{12}^{(n)} & = & y_2^{(n)} + hk_{21}/2 \\ k_{22}^{(n)} & = & -y_1^{(n)} - hk_{11}/2 + \sin(t_n + h/2) \\ k_{13}^{(n)} & = & y_2^{(n)} + hk_{22}/2 \\ k_{23}^{(n)} & = & -y_1^{(n)} - hk_{12}/2 + \sin(t_n + h/2) \\ k_{14}^{(n)} & = & y_2^{(n)} + hk_{23} \\ k_{24}^{(n)} & = & -y_1^{(n)} - hk_{13} + \sin(t_n + h). \end{eqnarray} $$

An example for a corresponding Octave implementation would be

function next = rhs(y, t)

    ## time step
    h = 2*pi/50;

    ## weights
    b = [1/6, 1/3, 1/3, 1/6];

    ## k_ij
    k = zeros(2, 4);
    k(1, 1) = y(2);
    k(2, 1) = (-y(1) + sin(t));
    k(1, 2) = y(2) + h * k(2, 1)/2;
    k(2, 2) = -y(1) - h * k(1, 1)/2 + sin(t + h/2);
    k(1, 3) = y(2) + h * k(2, 2)/2;
    k(2, 3) = -y(1) - h * k(1, 2)/2 + sin(t + h/2);
    k(1, 4) = y(2) + h * k(2, 3);
    k(2, 4) = -y(1) - h * k(1, 3) + sin(t + h);

    ## next time step
    next(1) = y(1) + h * b * k(1, :)';
    next(2) = y(2) + h * b * k(2, :)';

endfunction