# Rate of Convergence for Trapezoidal Method-System of Linear ODEs

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.

• What norm function are you using? Does it approximate an integral norm? The easiest to use would be the max norm, max(abs(u-uexact)). Or just consider the error at the end of the integration interval. – LutzL Sep 17 '15 at 22:36
• I was using the two norm. I went ahead and used the max norm and the error at the end of the integration interval and obtained the same rate of convergence as before. No luck so far. – Mario Sep 18 '15 at 0:40
• Order 2 makes one expect that for N=2^10=1024 the error is of the magnitude 1e-5 to 1e-6. However, your exact value does not have that accuracy. Try with more digits or use the formula as in the first component. – LutzL Sep 18 '15 at 8:36

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.