7 Implicit IVP Solvers
In some cases, IVPs can be difficult to solve because of the non-linearity of its terms, this is where Implicit Methods should be used to accommodate for these issues.
7.1 Backwards Euler Method
Consider the Euler method at the starting time \(t=t_0\). The value of the function \(y\) at \(t_1=t_0+h\) is approximated by \[ y(t_1) \approx Y_1=Y_0+hy'(t_0) \] and this gives an upper bound for a stable stepsize of \[ h_0=2\min\left(\frac{|\Re(\lambda_k)|}{|\lambda_k|^2}\right) \] in order to ensure that the Euler method is computationally stable. However, suppose that this modified slightly by using the gradient at \(y(t_1)\) rather than at \(y(t_0)\), in other words, suppose that the value of \(y\) at \(t_1\) is approximated by \[ y(t_1) \approx Y_1=Y_0+h\underline{\underline{y'(t_1)}}. \] This approach is known as the Backwards Euler Method and is an implicit procedure since the value of \(y'(t_1)\) is not known to begin with.
The general formulation is as follows: Consider the system of differential equations \[\boldsymbol{y}'=A\boldsymbol{y}+\boldsymbol{b}(t) \quad \text{with} \quad \boldsymbol{y}(t_0)=\boldsymbol{y}_0, \quad x \in [t_0,t_f].\] Discretise the interval \([t_0,t_f]\) into \(N\) equal subintervals, each with width \(h=\frac{t_f-t_0}{N}\). At the time step \(t=t_n=t_0+nh\), the backwards Euler method is \[\boldsymbol{Y}_{n+1}=\boldsymbol{Y}_n+h\boldsymbol{y}'(t_{n+1})=\boldsymbol{Y}_n+h\left[ A\boldsymbol{Y}_{n+1}+\boldsymbol{b}(t_{n+1}) \right].\] This can be rearranged to give \[(I-hA)\boldsymbol{Y}_{n+1}=\boldsymbol{Y}_n+h\boldsymbol{b}(t_{n+1}).\]
Rearranging further fives the basis for the Backwards Euler iteration which is \[\boldsymbol{Y}_{n+1}=(I-hA)^{-1}\left[\boldsymbol{Y}_n+h\boldsymbol{b}(t_{n+1})\right]\] whereas the standard Euler method in matrix form is \[\boldsymbol{Y}_{n+1}=(I+hA)^{-1}\boldsymbol{Y}_n+h\boldsymbol{g}(t_{n+1}).\] The Euler method requires explicit calculations using matrix multiplications but the backwards Euler method requires matrix inversion instead.
7.2 Stability of the Backwards Euler Method
Consider the initial value problem in its scalar form \[\frac{\mathrm{d} y}{\mathrm{d} t}=\lambda y+b(t) \quad \text{with} \quad y(0)=y_0.\] The backwards Euler method at the time \(t=t_{n+1}=t_0+(n+1)h\) gives \[Y_{n+1}=(1-h \lambda)^{-1}\left[Y_n+hg(t_{n+1})\right].\] This initial condition can be perturbed by adding a small parameter \(\varepsilon\neq 0\) to give the perturbed differential equation \[\frac{\mathrm{d} z}{\mathrm{d} t}=\lambda z+g(t) \quad \text{with} \quad z(0)=y_0+\varepsilon.\] The backwards Euler then yields \[Z_{n+1}=(1-h \lambda)^{-1}\left[Z_n+hg(t_{n+1})\right]\] The differential equations in \(Y\) and \(Z\) can be subtracted to give a perturbation term \(E\) where \[E_{n+1}=Z_{n+1}-Y_{n+1}=(1-h \lambda)^{-1}\left[Z_n-Y_n\right]=(1-h \lambda)^{-1}E_n.\] Notice that once again, the forcing function \(g(t)\) has been eliminated and therefore does not affect the stability of the backwards Euler method. The differential equation for \(E\) will have the initial condition \(E_0=Z_0-Y_0=\varepsilon\). This expression can be used to represent \(E_n\) in terms of \(\varepsilon\) recursively as: \[\begin{multline*} E_n=(1-h\lambda)^{-1}E_{n-1}=(1-h\lambda)^{-2}E_{n-2}\\ =\dots=(1-h\lambda)^{-(n-1)}E_1=(1-h\lambda)^{-n}E_0=(1-h\lambda)^{-n}\varepsilon. \end{multline*}\] \[\implies \quad E_n=(1-h\lambda)^{-n}\varepsilon.\] This means that the method is stable for stepsizes \(h\) that satisfy \(|1-h\lambda|>1\) and since \(\lambda<0\) for an asymptotically stable system, then this inequality is always satisfied. This means that the backwards Euler method is stable for all stepsizes \(h>0\), no matter how large.
Consider the differential equation \[ y'=-100y+100\sin(t) \quad \text{with} \quad y(0)=1. \] In this case, \(\lambda<0\) meaning that this differential equation is asymptotic stable. The maximum allowable stepsize for the Euler method is \(h_0=\frac{2}{|-100|}=0.02\). However, the backwards Euler method is stable for any stepsize \(h\) as seen below (very large stepsizes will still converge but they will not give any useful information).
The formulation presented above also holds for sets of differential equations in the same way with one difference. Instead of having \((1-h\lambda)^{-1}=\frac{1}{1-h\lambda}\), the procedure for systems will require the matrix inverse \((1-\lambda A)^{-1}\) or the MATLAB backslash operator can be used instead.
7.3 Order of Accuracy
The backwards Euler method is numerically stable for all values the stepsize \(h\) and has the same order of accuracy as the Euler method, i.e. the local truncation error is of \(\mathcal{O}\left(h^2\right)\) while the global integration error is of \(\mathcal{O}\left(h\right)\). However, this increased stability comes at a cost, the backwards Euler methods requires double the computational cost compared to the Euler method.
7.4 MATLAB Code
The following MATLAB code performs the Backwards Euler iteration for the system \(\boldsymbol{y}'=A\boldsymbol{y}+\boldsymbol{b}(t)\) subject to \(\boldsymbol{y}(0)=\boldsymbol{y}_0\) where \[A=\begin{pmatrix} -7 & -2 & 1 \\ 2 & -1 & -9 \\ 0 & 0 & -5 \end{pmatrix}, \quad \boldsymbol{b}(t)=\begin{pmatrix} \sin(t) \\ 0 \\ 2 \end{pmatrix}, \quad \boldsymbol{y}_0=\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}.\]
function IVP_Back_Euler
%% Solve a set of first order IVPs using Backwards Euler
% This code solves a set of IVP when written in the form
% dydt=A*y+b(t) on the interval [t0,tf] subject to the
% initial conditions y(0)=y0. The output will be the graph
% of the solution(s) and the vector value at the final
% point tf.
%% Lines to change:
% Line 25 : t0 - Start time
% Line 28 : tf - End time
% Line 31 : N - Number of subdivisions
% Line 34 : A - Matrix A
% Line 37 : b - Forcing vector b(t)
% Line 40 : y0 - Vector of initial values
% Line 106+ : Which functions to plot, remembering to assign
% a colour, texture and legend label
%% Set up input values
% Start time
t0=0;
% End time
tf=1;
% Number of subdivisions
N=5000;
% Matrix A
A=[-7,-2,1;2,-1,-9;0,0,-5];
% Vector b, which can be a function of t in general
b=@(t) [sin(t);0;2];
% Column vector initial values y0=y(t0)
y0=[0;1;0];
%% Set up IVP solver parameters
% T = Vector of times t0,t1,...,tN.
% This is generated using linspace which splits the
% interval [t0,tf] into N+1 points (or N subintervals)
T=linspace(t0,tf,N+1);
% Stepsize
h=(tf-t0)/N;
% Number of differential equations
K=length(y0);
%% Perform the Euler iteration
% Y = Solution matrix
% The matrix Y will contain K rows and N+1 columns. Every
% row corresponds to a different IVP and every column
% corresponds to a different time. So the matrix Y will
% take the following form:
% y_1(t_0) y_1(t_1) y_1(t_2) ... y_1(t_N)
% y_2(t_0) y_2(t_1) y_2(t_2) ... y_2(t_N)
% ...
% y_K(t_0) y_K(t_1) y_K(t_2) ... y_K(t_N)
Y=zeros(K,N+1);
% The first column of the vector Y is the initial vector y0
Y(:,1)=y0;
% Set the current time t to be the starting time t0 and the
% current value of the vector y to be the strtaing values y0
y=y0;
for n=2:1:N+1
t=T(n); % Update the new time
y=(eye(K,K)-h*A)\(y+h*b(t));% Find y at the current step
Y(:,n)=y; % Replace row n in Y with y
end
%% Setting plot parameters
% Clear figure
clf
% Hold so more than one line can be drawn
hold on
% Turn on grid
grid on
% Setting font size and style
set(gca,'FontSize',20,'FontName','Times')
set(legend,'Interpreter','Latex')
% Label the axes
xlabel('$t$','Interpreter','Latex')
ylabel('$\mathbf{y}(t)$','Interpreter','Latex')
% Plot the desried solutions. If all the solutions are
% needed, then consider using a for loop in that case
plot(T,Y(1,:),'-b','LineWidth',2,'DisplayName','$y_1(t)$')
plot(T,Y(2,:),'-r','LineWidth',2,'DisplayName','$y_1(t)$')
plot(T,Y(3,:),'-k','LineWidth',2,'DisplayName','$y_1(t)$')
end7.5 Stiff Differential Equations
Stiff sets of differential equations with a large value of the total computational cost \(N_0\) can be very difficult to solve numerically using explicit methods but implicit methods can work very well. MATLAB has its very own built-in stiff differential equation solver under the command ode15s and this can be implemented exactly as ode45. This solves sets of differential equations implicitly using numerical differentiation of orders 1 to 5.
Consider the set of differential equations on the interval \([0,3500]\) \[\begin{align*} & \frac{\mathrm{d} y_1}{\mathrm{d} t}=y_2 & y_1(0)=2 \\ & \frac{\mathrm{d} y_2}{\mathrm{d} t}=1000(1-y_1^2)y_2-y_1 & y_2(0)=0. \end{align*}\]
This is a very stiff set of differential equations, solving this using ode45 takes upwards of 92 seconds while solving using the stiff solver ode15s requires a mere 0.233 seconds (depending on you machine). The result of solving this differential equation is shown below for \(y_1(t)\) only since \(y_2(t)\) takes very large values and this distorts the graphical interpretation.
Using the stiff solver optimises the stepsizes for stiff regions. Particularly, if a region is deemed to be considerably “stiff”, the ode15s will use smaller stepsizes to solve the problem but if there is a region where the differential is not “stiff”, then larger stepsizes will be used. Therefore, ode15s usually requires fewer grid points overall, for instance to solve the above set of differential equations, ode15s only requires 1,836 grid points while ode45 requires 7,820,485 grid points, that is over 4,200 times more grid points than ode15s. This just goes to show that stiff differential need implicit methods, even though the cost for every step is greater than that of an explicit method, fewer steps are required in total.
An alternative stiff differential equation solver is ode23s which achieves that same outcome as ode15s but with a lower accuracy and more grid points using only second and third order methods.