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.