[Math] Rate of Convergence for Trapezoidal Method-System of Linear ODEs

numerical methodsordinary differential equations

I have been given a system of two linear ODEs that I am to solve using various analytical and numerical methods:

$u_{1}' = u_{1} $
$u_{2}' = u_{1} – u_{2}$

I have solved for the solution of these equations using analytical methods (Duhamel's principle, exponential matrix, etc.) and was able to solve numerical using Euler's method and was able to show that the rate of convergence is one.

One of the methods I am supposed to use is the trapezoidal method and I need to show numerically the rate of convergence of two. Using the trapezoidal formula, I rewrite the equations to solve for the next step. Since both ODEs are linear, I can explicitly solve for the subsequent step.

$u_{1}^{n+1} = u_{1}^{n}*\frac{(1+0.5*dt)}{(1-0.5*dt)}$

$u_{2}^{n+1} = \frac{(u_{2}^{n} + 0.5*dt*(u_{1}^{n}-u_{2}^{n}+u_{1}^{n+1})}{( 1+0.5*dt)}$

To determine the rate of convergence, I wrote a MATLAB script to calculate the slope of the logarithm of the error over the logarithm of the number of interval. 

uexact = [exp(1),1.5431];

N = 2.^[5:10];
T = 1.0;

for i = 1:length(N)

u1(1) = 1;
u2(1) = 1;
dt = T/N(i);

for j = 1:N(i)

    u1(j+1) = u1(j) * (1+0.5*dt)/(1-0.5*dt);
    u2(j+1) = (u2(j) + 0.5*dt*(u1(j)-u2(j)+u1(j+1))) / ( 1+0.5*dt);

end

u = [u1(end),u2(end)];

err(i)= norm(u-uexact);

clear u1 u2


end

loglog(N,err)


m = -(log(err(end))-log(err(1)))...
/(log(N(end))-log(N(1))

However, when I do this, I calculate a rate of convergence of 0.7204. I have no clue what I am doing wrong. Any obvious mistakes or mistakes in the thinking process?

I went ahead for the case where I only have one ODE and used the trapezoidal method and calculated a rate of convergence of 2. For the system however, I've had no luck.

Thank you for the help.

Best Answer

Your exact value is in the second component not exact enough. With N ranging from 32 to 1024, one would expect for an order 2 method and a tame ODE system errors from 1e-3 to 1e-6. However, if the supposedly exact value has an error of magnitude 1e-4 to 1e-5, then the integration error at the last iterations get unnecessarily perturbed.

\begin{align} u_1(t) &= \exp(t)\\ (u_2(t)·\exp(t))' &= u_1(t)·\exp(t) = \exp(2·t)\\ u_2(t)·\exp(t) &= C + 0.5·\exp(2·t) \implies C = 0.5\\ u_2(t) &= (\exp(t)+\exp(-t))/2 = \cosh(t) \end{align}

Related Question