Implementing the $\theta$-method for solving an IVP

Question

I also posted this question on the Math SE site, but it was suggested I post here also.

I want to implement the $\theta$-method to solve an IVP in MATLAB. The $\theta$-method is:

$y_{j+1} = y_j + h[\theta f(t_j, y_j) + (1 - \theta)f(t_{j+1}, y_{j+1})]$ for $\theta \in [0, 1]$.

I want to use it to solve an IVP of this form:

$y_1\prime = f_1(y_1, y_2) \;\;\;\;\; y_2\prime = f_2(y_1, y_2)$

$t \in [a, b]$

$y_1(a) = \alpha \;\;\;\;\; y_2(a) = \alpha \;\;\;\;\;\;$ (The initial conditions)

I'm not sure how I would implement this, since I'm used to implementing methods in which there is only one DE, such as $y\prime=f(t, y)$ (like Euler's method or the Runge-Kutta).

Here is what I did so far in MATLAB:

function thetaMeth(a, b, h, f, y)

    j = 1;
    theta = 0.5
    for i = a:h:b
        y(j+1) = y(j) + h*(theta*f(j, y(j)) + (1 - theta)*f(j+1, y(j+1)));
        fprintf('t = %f; y(%d) = %f\n', i, j, y(j));
    end
        j = j + 1;
end

My parameter f would be for the $y_1\prime$ and $y_2\prime$ equations (it could possibly be taken in as a cell array of functions). I'm not sure what to do with the y parameter, because the solution equations $y_1$ and $y_2$ are not given (I'm not really sure if I should even have y as a parameter). This is an implicit method, and you need to have $y_{j+1}$ calculated beforehand somehow, because it appears in the LHS and RHS.

Anyone have any pointers as to how I can implement this? It doesn't have to be in MATLAB, I just need an algorithm that can be translated into code, since I'm not sure how to implement the algorithm myself.

Answer 1

First of all, the formula that you use for the $\theta$-method is correct, but beware that you mix the indices for the time advancing with the indices of the components. Write for example $y^{(i)}_j$ for the component $j$ of the solution at time $t_i$. This can avoid confusions.

If the system of ODEs that you want to solve is

$$ \mathbf{y}\prime = f(t,\mathbf{y}) $$

i.e. written in components

$$ \left\{\begin{array}{c} y_1\prime = f_1(t,\mathbf{y}) \\ y_2\prime = f_2(t,\mathbf{y}) \end{array} \right.$$

with $\mathbf{y}= \left( \begin{array}{c} y_1 \\ y_2 \end{array} \right)$ and $f(t,y) = \left( \begin{array}{c} f_1(t,\mathbf{y}) \\ f_2(t,\mathbf{y}) \end{array} \right)$. As you mention, the scheme is $$\mathbf{y}^{(j+1)} = \mathbf{y}^{(j)} + h \left( \theta f(t_j, \mathbf{y}^{(j)}) + (1-\theta)f(t_{j+1}, \mathbf{y}^{(j+1)})\right)$$

which is implicit. The fact that it is implicit scheme implies that (in general) you cannot find directly $\mathbf{y}^{(j+1)}$ from $\mathbf{y}^{(j)}$, but you need to solve non-linear equations. Indeed, you can rewrite the scheme as $$G(\mathbf{y}^{(j+1)}) =\mathbf{y}^{(j+1)} - \mathbf{y}^{(j)} - h \left( \theta f(t_j, \mathbf{y}^{(j)}) + (1-\theta)f(t_{j+1}, \mathbf{y}^{(j+1)})\right) = \mathbf{0}$$.

So, at each time step of your algorithm, you need to solve this non-linear equation. This can be done for example with the Matlab function fzero (applied to each component of the equation). That is what you should do in your for loop.

For the implementation, you can use an anonymous function for f, it supports vectorial valued functions. Instead of y, take alpha as argument. Then, y is the output of your function (possible with the times to which the approximations correspond).

Hope it helps!