Section 19.5 Solving systems of differential equations
Conceptually, any of the above methods can be directly applied to systems of \(m\) differential equations, just by interpreting \(y\) as a vector function \(\mathbf y\text{.}\) The notation becomes awkward: now \(y_2\) refer to the second component of \(\mathbf y\text{,}\) while the value \(\mathbf y(t_2)\) might be denoted \(\mathbf y_2\text{.}\) There are also a few changes in Matlab notation regarding y
: it will now be a matrix, where rows correspond to components and columns correspond to points of evaluation. The right-hand side function f
should return a column vector of the same dimesion \(m\text{.}\)
Example 19.5.1. Implement trapezoidal method for ODE system.
Solve the system
with initial condition \(y(0) = \begin{pmatrix}1 \\ 0\end{pmatrix}\text{,}\) on the interval \([0, 10]\) using the trapezoidal method with the step \(h=0.01\text{.}\) Plot the solution.
f = @(t, y) [t-3*y(2); 2*y(1)]; h = 0.01; t = 0:h:10; y0 = [1; 0]; y = y0*ones(size(t)); for k = 1:numel(t)-1 pred = y(:, k) + h*f(t(k), y(:, k)); y(:, k+1) = y(:, k) + h/2*(f(t(k), y(:, k)) + f(t(k+1), pred)); end plot(t, y)
The most visible difference is the notation y(:, k)
which takes the entire \(k\)th column of the matrix y
. Its mathematical meaning is the approximation to \(\mathbf y(t_k)\text{.}\) The initial condition y0
is a column vector, which makes the computation y = y0*ones(size(t));
an outer product: it creates a matrix in which every column is equal to y0
. This is convenient because we need the first column to be equal, and the other columns will be rewritten anyway.
The plot displayed by Example 19.5.1 is a time series plot in which one axis represents the independent variable (usually time) and the other axis shows the components of \(\mathbf y\text{.}\) The other option we have is to make a phase plot in which each component of \(\mathbf y\) gets its own axis, and there is no time axis. This is achieved with
plot(y(1,:), y(2,:))
The phase plot may be preferable when the system does not involve time directly (an autonomous system). Try both kinds of plots for the autonomous system
f = @(t, y) [y(2)-y(2)^3; 2*y(1)];