# Numerical estimates for the convergence order of trapezoidal-like Runge-Kutta methods

I have to calculate approximations of the solution with the method $$y^{n+1}=y^n+h \cdot [\rho \cdot f(t^n,y^n)+(1-\rho) \cdot f(t^{n+1},y^{n+1})] ,\quad n=0,\ldots,N-1 \\ y^0=y_0$$ for various values of $\rho$ and the errors for uniform partitions with $N=64, 128, \ldots, 4096,8192$ subintervals. Determine the value of the parameter $\rho \in [0,1)$ such that the method has maximum order, let $\rho=\rho_m$.

PS: The order accuracy can be computed numerically as follows.

Let $E(N_1)$ (respectively $E(N_2)$) the error of the numerical method for $N_1$ ( respectively $N_2$) subintervals, and let's suppose that $E(N_k) \approx Ch_k^p$ where the constant $C$ is independent of $h_k$ and $N_k, k=1,2$.Then $\frac{E(N_1)}{E(N_2)} \approx \frac{Ch_1^p}{Ch_2^p} \Rightarrow p \approx \frac{\log\frac{E(N_1)}{E(N_2)} }{\log\frac{h_1}{h_2}}$.

How can we find the parameter $\rho$?

EDIT

The ode is $$y'=y+4\pi \cos (4\pi t)y, t\in [0,1]\\ y(0)=1$$

The exact solution is $$y(t)=e^{t+\sin (4\pi t)}$$

EDIT 2:

I implemented the method and I got the following results:

     rho= 0.000000
p1(2.000000)= 1.083617
p1(3.000000)= 1.004083
p1(4.000000)= 1.000023
p1(5.000000)= 0.999565
p1(6.000000)= 0.999678
p1(7.000000)= 0.999814
p1(8.000000)= 0.999901
p1_max(1.000000) = 1.083617
rho= 0.100000
p1(2.000000)= 1.058614
p1(3.000000)= 0.977502
p1(4.000000)= 0.987617
p1(5.000000)= 0.993574
p1(6.000000)= 0.996735
p1(7.000000)= 0.998355
p1(8.000000)= 0.999175
p1_max(2.000000) = 1.058614
rho= 0.200000
p1(2.000000)= 1.032348
p1(3.000000)= 0.951871
p1(4.000000)= 0.975554
p1(5.000000)= 0.987725
p1(6.000000)= 0.993855
p1(7.000000)= 0.996926
p1(8.000000)= 0.998463
p1_max(3.000000) = 1.032348
rho= 0.300000
p1(2.000000)= 1.034457
p1(3.000000)= 0.927331
p1(4.000000)= 0.963963
p1(5.000000)= 0.982096
p1(6.000000)= 0.991082
p1(7.000000)= 0.995550
p1(8.000000)= 0.997777
p1_max(4.000000) = 1.034457
rho= 0.400000
p1(2.000000)= 1.017292
p1(3.000000)= 0.938286
p1(4.000000)= 0.953174
p1(5.000000)= 0.976886
p1(6.000000)= 0.988522
p1(7.000000)= 0.994282
p1(8.000000)= 0.997146
p1_max(5.000000) = 1.017292
rho= 0.500000
p1(2.000000)= 0.992201
p1(3.000000)= 0.999250
p1(4.000000)= 0.989454
p1(5.000000)= 0.972753
p1(6.000000)= 0.986534
p1(7.000000)= 0.993307
p1(8.000000)= 0.996663
p1_max(6.000000) = 0.999250
rho= 0.600000
p1(2.000000)= 0.977841
p1(3.000000)= 0.988846
p1(4.000000)= 0.995149
p1(5.000000)= 0.997443
p1(6.000000)= 0.998773
p1(7.000000)= 0.999399
p1(8.000000)= 0.999699
p1_max(7.000000) = 0.999699
rho= 0.700000
p1(2.000000)= 0.953567
p1(3.000000)= 0.978911
p1(4.000000)= 0.989106
p1(5.000000)= 0.994622
p1(6.000000)= 0.997315
p1(7.000000)= 0.998651
p1(8.000000)= 0.999326
p1_max(8.000000) = 0.999326
rho= 0.800000
p1(2.000000)= 0.934816
p1(3.000000)= 0.966643
p1(4.000000)= 0.982599
p1(5.000000)= 0.991290
p1(6.000000)= 0.995600
p1(7.000000)= 0.997797
p1(8.000000)= 0.998897
p1_max(9.000000) = 0.998897
rho= 0.900000
p1(2.000000)= 0.910235
p1(3.000000)= 0.953068
p1(4.000000)= 0.975759
p1(5.000000)= 0.987650
p1(6.000000)= 0.993752
p1(7.000000)= 0.996861
p1(8.000000)= 0.998430
p1_max(10.000000) = 0.998430
rho= 1.000000
p1(2.000000)= 0.888528
p1(3.000000)= 0.939543
p1(4.000000)= 0.968448
p1(5.000000)= 0.983873
p1(6.000000)= 0.991846
p1(7.000000)= 0.995897
p1(8.000000)= 0.997945
p1_max(11.000000) = 0.997945

rho=0.000000, p_max=1.083617


p1() is given from the relation:

p1(i)=(log(E(N(i-1))/E(N(i)))/(log((1./N(i-1))/(1./N(i)))));

p1_max() is the maximum p1() for a specific rho.

p_max is the maximum values of all p1_max(). We get this values for rho=0.000000

But this cannot be trues, because that would mean that backward Euler method has the maximum order.

EDIT 3: I get the following $y_{N(i)}$ and errors:

     rho= 0.000000
i=1.000000 N(1.000000)=1024.000000 y(1) = 2.789277 , error=0.144219

i=2.000000 N(2.000000)=2048.000000 y(1) = 2.753557 , error=0.071584

i=3.000000 N(3.000000)=4096.000000 y(1) = 2.735865 , error=0.035663

i=4.000000 N(4.000000)=8192.000000 y(1) = 2.727060 , error=0.017799

rho= 0.250000
i=1.000000 N(1.000000)=1024.000000 y(1) = 2.735357 , error=0.089452

i=2.000000 N(2.000000)=2048.000000 y(1) = 2.726813 , error=0.044631

i=3.000000 N(3.000000)=4096.000000 y(1) = 2.722546 , error=0.022292

i=4.000000 N(4.000000)=8192.000000 y(1) = 2.720413 , error=0.011140

rho= 0.500000
i=1.000000 N(1.000000)=1024.000000 y(1) = 2.682558 , error=0.035826

i=2.000000 N(2.000000)=2048.000000 y(1) = 2.700348 , error=0.017960

i=3.000000 N(3.000000)=4096.000000 y(1) = 2.709297 , error=0.008992

i=4.000000 N(4.000000)=8192.000000 y(1) = 2.713785 , error=0.004499

rho= 0.750000
i=1.000000 N(1.000000)=1024.000000 y(1) = 2.630851 , error=-0.016687

i=2.000000 N(2.000000)=2048.000000 y(1) = 2.674158 , error=-0.008433

i=3.000000 N(3.000000)=4096.000000 y(1) = 2.696117 , error=-0.004239

i=4.000000 N(4.000000)=8192.000000 y(1) = 2.707173 , error=-0.002125

rho= 1.000000
i=1.000000 N(1.000000)=1024.000000 y(1) = 2.580212 , error=-0.068114

i=2.000000 N(2.000000)=2048.000000 y(1) = 2.648241 , error=-0.034552

i=3.000000 N(3.000000)=4096.000000 y(1) = 2.683006 , error=-0.017401

i=4.000000 N(4.000000)=8192.000000 y(1) = 2.700579 , error=-0.008732


EDIT 4: I get the following results:

rho= 0.000000
i=1.000000 N(1.000000)=1024.000000000  y(1) =    2.826726525  , error=   0.144219272

i=2.000000 N(2.000000)=2048.000000000  y(1) =    2.771919117  , error=   0.071584132

i=3.000000 N(3.000000)=4096.000000000  y(1) =    2.744956198  , error=   0.035662690

i=4.000000 N(4.000000)=8192.000000000  y(1) =    2.731583189  , error=   0.017799247

rho 0.250000
i=1.000000 N(1.000000)=1024.000000000  y(1) =    2.771959612  , error=   0.089452359

i=2.000000 N(2.000000)=2048.000000000  y(1) =    2.744966075  , error=   0.044631089

i=3.000000 N(3.000000)=4096.000000000  y(1) =    2.731585628  , error=   0.022292121

i=4.000000 N(4.000000)=8192.000000000  y(1) =    2.724924189  , error=   0.011140247

rho= 0.500000
i=1.000000 N(1.000000)=1024.000000000  y(1) =    2.718333217  , error=   0.035825965

i=2.000000 N(2.000000)=2048.000000000  y(1) =    2.718294675  , error=   0.017959690

i=3.000000 N(3.000000)=4096.000000000  y(1) =    2.718285040  , error=   0.008991533

i=4.000000 N(4.000000)=8192.000000000  y(1) =    2.718282631  , error=   0.004498690

rho= 0.750000
i=1.000000 N(1.000000)=1024.000000000  y(1) =    2.665819881  , error=  -0.016687372

i=2.000000 N(2.000000)=2048.000000000  y(1) =    2.691901512  , error=  -0.008433474

i=3.000000 N(3.000000)=4096.000000000  y(1) =    2.705054009  , error=  -0.004239498

i=4.000000 N(4.000000)=8192.000000000  y(1) =    2.711658463  , error=  -0.002125479

rho= 1.000000
i=1.000000 N(1.000000)=1024.000000000  y(1) =    2.614392911  , error=  -0.068114342

i=2.000000 N(2.000000)=2048.000000000  y(1) =    2.665783224  , error=  -0.034551762

i=3.000000 N(3.000000)=4096.000000000  y(1) =    2.691892114  , error=  -0.017401394

i=4.000000 N(4.000000)=8192.000000000  y(1) =    2.705051630  , error=  -0.008732312

p1(2.000000)= 0.991846
p1(3.000000)= 0.995897
p1(4.000000)= 0.997945
rho=0.000000, p_max=1.010552


when I execute the following: http://pastebin.com/TyYQpMA3

and the method is this: http://pastebin.com/pxTxE2Rf

Have I maybe done something wrong at the calculation of the error E?

I find the exact solution here: http://pastebin.com/7HX65U8W

EDIT 5: The results for the changed version of the function of the exact solution are the following:

     rho= 0.000000
i=1.000000 N(1.000000)=1024.000000000  y(1) =  2.82672652466  , error=0.1084446962030

i=2.000000 N(2.000000)=2048.000000000  y(1) =  2.77191911745  , error=0.0536372889879

i=3.000000 N(3.000000)=4096.000000000  y(1) =  2.74495619768  , error=0.0266743692175

i=4.000000 N(4.000000)=8192.000000000  y(1) =  2.73158318927  , error=0.0133013608113

rho= 0.250000
i=1.000000 N(1.000000)=1024.000000000  y(1) =  2.77195961188  , error=0.0536777834160

i=2.000000 N(2.000000)=2048.000000000  y(1) =  2.74496607477  , error=0.0266842463094

i=3.000000 N(3.000000)=4096.000000000  y(1) =  2.73158562815  , error=0.0133037996933

i=4.000000 N(4.000000)=8192.000000000  y(1) =  2.72492418927  , error=0.0066423608061

rho= 0.500000
i=1.000000 N(1.000000)=1024.000000000  y(1) =  2.71833321749  , error=0.0000513890326

i=2.000000 N(2.000000)=2048.000000000  y(1) =  2.71829467535  , error=0.0000128468916

i=3.000000 N(3.000000)=4096.000000000  y(1) =  2.71828504016  , error=0.0000032117000

i=4.000000 N(4.000000)=8192.000000000  y(1) =  2.71828263138  , error=0.0000008029236

rho= 0.750000
i=1.000000 N(1.000000)=1024.000000000  y(1) =  2.66581988059  , error=-0.0524619478724

i=2.000000 N(2.000000)=2048.000000000  y(1) =  2.69190151155  , error=-0.0263803169049

i=3.000000 N(3.000000)=4096.000000000  y(1) =  2.70505400927  , error=-0.0132278191846

i=4.000000 N(4.000000)=8192.000000000  y(1) =  2.71165846267  , error=-0.0066233657940

rho= 1.000000
i=1.000000 N(1.000000)=1024.000000000  y(1) =  2.61439291112  , error=-0.1038889173386

i=2.000000 N(2.000000)=2048.000000000  y(1) =  2.66578322391  , error=-0.0524986045516

i=3.000000 N(3.000000)=4096.000000000  y(1) =  2.69189211409  , error=-0.0263897143722

i=4.000000 N(4.000000)=8192.000000000  y(1) =  2.70505163034  , error=-0.0132301981164

p1(2.000000)= 0.984691
p1(3.000000)= 0.992303
p1(4.000000)= 0.996141
rho=0.500000, p_max=2.000000


EDIT 5: That is the graph that I get for the error in respect to $N$:

• does the solution $y$ satisfy some ODE, possibly of the form $y'=f(t,y)$? Commented Mar 30, 2015 at 23:53
• @LutzL I added the differential equation, the interval and the exact solution. To implement this method do we not have to take a specific $\rho$? Commented Mar 31, 2015 at 11:31
• Yes, you have to discretize the domain of $ρ$ too. I'd propose $ρ_k=k/10$, $k=0,1,...,10$. The solution of the implicit scheme is easy to directly encode for linear ODE. Commented Mar 31, 2015 at 11:44
• In my original reading I thought you were asking about how t find $p$. $\rho$ is in the definition of your solution technique: it is the mix of $f(t^n,y^n)$ and $f(t^{n+1},y^{n+1})$ that you use. You choose it. $\rho=1$ is the forward Euler method, $rho=0$ is the backward Euler method, and values in between are methods in between. $p$ is the order of accuracy of the resulting solution. It is defined by what was the last equation. Commented Mar 31, 2015 at 14:12
• I can't see anything wrong, the Newton iteration should be correct in the first step and finish after the second, independent of TOL, since the equation is linear. The only thing you could try is to iterate for k in y(n+1)=y(n)+h*k with $k=\rho f(t_n,y_n)+(1-\rho) f(t_{n+1},y_n+hk)$, this might avoid some rounding or truncation errors. (And get rid of the condition if (k>Kmax || dx < TOL) after the loop.) Commented Apr 1, 2015 at 7:52

Building upon the answer of Ross Millikan: If the exact solution is not known, then one can only use it indirectly. The exact solution is $y(T)$ and the approximative solutions give results $$Y(N)=y(T)+Ch_N^p+O(h_N^{p+1})$$ Comparing the solution with $N$ steps of step-size $h_N=T/N$ and $Y(N)=y_N$ the last iterate, with the solution with $N/2$ steps gives $$Y(N/2)-Y(N)=Ch_N^p(2^p-1)+O(h_N^{p+1})$$ Doubling the number of steps gives $$\frac{Y(N/2)-Y(N)}{Y(N)-Y(2N)}=\frac{C+O(h_N)}{C+O(2^ph_N)}·2^p\approx 2^p$$ so that again the logarithm serves to get an approximate value for $p$ $$p\approx \frac{\log\bigl(Y(N/2))-Y(N)\bigr)-\log\bigl((N)-Y(2N)\bigr)}{\log(2)}.$$

## Numerical experiments

Python script

from math import pi, cos, log

def c(t):
return 1+4*pi*cos(4*pi*t);

def integrate(N, rho):
h=1.0/N;
y=1.0
for k in range(N):
t=k*h;
y = y * (1+rho*h*c(t))/(1-(1-rho)*h*c(t+h));
return y;

for k in range(5):
rho = k/4.0;
print "    rho = ", rho;
N=1024;
ylast = 0
for i in range(4):
y = integrate(N, rho);
if(i>0):
err = ylast-y
print "    i=", i, " N=",N,"\ty(1)=",y,"\terr(",i-1,")=",err,"\tlog2|err|=",log(abs(err))/log(2);
else:
print "    i=", i, " N=",N,"\ty(1)=",y
ylast = y;
N *= 2;
print "    --";


Cursory results

rho =  0.0
i= 0  N= 1024       y(1)= 2.82672652466
i= 1  N= 2048       y(1)= 2.77191911745     err( 0 )= 0.0548074072152       log2|err|= -4.18948530333
i= 2  N= 4096       y(1)= 2.74495619768     err( 1 )= 0.0269629197704       log2|err|= -5.21287945772
i= 3  N= 8192       y(1)= 2.73158318927     err( 2 )= 0.0133730084062       log2|err|= -6.22453213762
--
rho =  0.25
i= 0  N= 1024       y(1)= 2.77195961188
i= 1  N= 2048       y(1)= 2.74496607477     err( 0 )= 0.0269935371066       log2|err|= -5.21124215659
i= 2  N= 4096       y(1)= 2.73158562815     err( 1 )= 0.0133804466161       log2|err|= -6.22372991834
i= 3  N= 8192       y(1)= 2.72492418927     err( 2 )= 0.0066614388872       log2|err|= -7.22995044802
--
rho =  0.5
i= 0  N= 1024       y(1)= 2.71833321749
i= 1  N= 2048       y(1)= 2.71829467535     err( 0 )= 3.85421409912e-05     log2|err|= -14.6632037598
i= 2  N= 4096       y(1)= 2.71828504016     err( 1 )= 9.6351915686e-06      log2|err|= -16.6632552186
i= 3  N= 8192       y(1)= 2.71828263138     err( 2 )= 2.40877642854e-06     log2|err|= -18.6632680738
--
rho =  0.75
i= 0  N= 1024       y(1)= 2.66581988059
i= 1  N= 2048       y(1)= 2.69190151155     err( 0 )= -0.0260816309674      log2|err|= -5.26082210107
i= 2  N= 4096       y(1)= 2.70505400927     err( 1 )= -0.0131524977203      log2|err|= -6.24851938976
i= 3  N= 8192       y(1)= 2.71165846266     err( 2 )= -0.00660445339054     log2|err|= -7.24234512113
--
rho =  1.0
i= 0  N= 1024       y(1)= 2.61439291112
i= 1  N= 2048       y(1)= 2.66578322391     err( 0 )= -0.051390312787       log2|err|= -4.28235975665
i= 2  N= 4096       y(1)= 2.69189211409     err( 1 )= -0.0261088901794      log2|err|= -5.25931505601
i= 3  N= 8192       y(1)= 2.70505163034     err( 2 )= -0.0131595162558      log2|err|= -6.24774973324
--


The errors show the halving behavior except for $\rho=0.5$, where the error gets quartered. The dyadic logarithm of the error also shows that behavior by the arithmetic progressions by about $-p=-1$ except in the middle where it progresses by $-p=-2$. With a denser set of $\rho$ values one sees the error getting progressively smaller from both sides towards $\rho=0.5$ in the middle without losing the order $p=1$, with then a jump to much smaller errors at $\rho=0.5$ with order $p=2$.

• I made some changes and now I get the same $y(1)$ as you got. Is err(i) the local error? Also what results do you get for the order $$p \approx \frac{\log \left (\frac{E(N_1)}{E(N_2)}\right )}{\log \left (\frac{h_1}{h_2}\right )}$$ ? Commented Apr 1, 2015 at 22:14
• err(i) is the difference between the last two numerical solutions. Approximately $C·(h/2)^p-C·h^p=D·h^p$. The log2|err(i)| of that has thus, with $h=1/N=2^{-n}$, the approximate form $\log_2|D|+p·\log_2(h)=\log_2|D|-p·n$. $-p$ is thus the step of the arithmetic progression of log2|err(i)|. In the end this is the formula I gave above that does not use the exact value. However, with the exact value the only change is that the constant is $\log_2|C|$ instead of $\log_2|D|$. Commented Apr 2, 2015 at 7:30
• I edited my post... Could you take a look at it? Commented Apr 2, 2015 at 8:32
• @evinda: Since the solutions are the same, and for $\rho=0.5$ get very close to the correct result $\exp(1)$ and the difference between discretizations is as expected.... there has to be something wrong with your exact solution. And indeed, you compute as $y_{exact}(t)$ always the value $y(t-h)$ by evaluating y_ex(i+1) at t(i). This of course has automatically an error $O(h)$ that is impossible to reduce by the numerical method. Commented Apr 2, 2015 at 9:02
• Yes, this is looking good. You could invest 20 min, but not more, into finding out if it is possible to label the scales in powers of 2. Commented Apr 2, 2015 at 18:40

Your last equation is exactly what you want to find $p$. If you have an analytic solution and have done numeric integration with $N_1,N_2$ steps, you can find the error for each integration. You also know the stepsizes $h_1,h_2$, so you have the whole right hand side.

• To do numeric integration with $N_1,N_2$ steps do we not have to take a specific $\rho$ ? Commented Mar 31, 2015 at 11:23
• Yes, you do. I think the point of the problem is to choose a $\rho$ in the range $[0,1]$, do the integration with varying numbers of steps (at the same $\rho$) and see if the accuracy $p$ improves. Commented Mar 31, 2015 at 14:15
• math.stackexchange.com/questions/1197371/… Could you take a look at my question? Commented Apr 30, 2015 at 17:08

So as an addendum to my comment I assume $y'=f(t,y)$ and also I assume $f$ is lipschitz continuous so there exists a unique solution, now by taylor expansion we get that $$y(t^{n+1}) = y(t^n)+hf(t^n,y(t^n))+(h^2/2)y''(\xi_1)$$ and $$y(t^{n+1}) = y(t^n)+hf(t^{n+1},y(t^{n+1}))-(h^2/2)y''(\xi_2)$$

and note that we can write $y(t^{n+1})=\rho y(t^{n+1})+(1-\rho)y(t^{n+1})$, subtracting this from our method we get: $$|y^{n+1}-y(t^{n+1})|=|y^{n+1}-(\rho y(t^{n+1})+(1-\rho)y(t^{n+1}))|$$ $$\le|y^n-y(t^n)|+h\rho|f(t^n,y^n)-f(t^n,y(t^n))|+h(1-\rho)|f(t^{n+1},y^{n+1})-f(t^{n+1},y(t^{n+1}))|$$ $$+h^2\|y''\|_\infty$$ $$\le(1+h\rho L)|y^n-y(t^n)|+h(1-\rho)L|y^{n+1}-y(t^{n+1})|+h^2\|y''\|_\infty$$ so $$|y^{n+1}-y(t^{n+1})|\le \frac{1+h\rho L}{1-hL+\rho hL}|y^n-y(t^n)|+\frac{h^2}{1-hL+\rho hL}\|y''\|_\infty,$$ iterating and letting $\gamma = \frac{1+h\rho L}{1-hL+h\rho L}$, and using the fact that $1+x+...+x^n=\frac{x^{n+1}-1}{x-1}$ we get $$|y^{n+1}-y(t^{n+1})|\le\frac{\gamma^{n+1}-1}{\gamma -1}\frac{h^2}{1-hL+h\rho L}\|y''\|_\infty$$

• @LutzL cheers edited :) Commented Mar 31, 2015 at 7:27
• @LutzL also $\gamma = O(1)$, so $1+\gamma+...+\gamma^n=\frac{\gamma^{n+1}-1}{\gamma-1}=O(1)$ also, so it doesn't affect the order Commented Mar 31, 2015 at 7:30
• That is true, but strange. For $ρ=1$ the method is Euler forward, for $ρ=0$ Euler backward, both of global error order $1$. It would be truly revolutionary if the $O(h^2)$ bound were valid. Commented Mar 31, 2015 at 7:36
• My earlier remark remains: Since $γ=1+Lh+O(h^2)$ one gets $γ^n=exp(n*(Lh+O(h^2))=exp(Lt_n+O(t_nh))$ so that the numerator is in the first order independent of $h$ and the denominator $γ-1=Lh+O(h^2)$ is in the first order proportional to $h$. Your bound is thus $O(e^{LT}h)$ for any fixed end time $T\gg h$. Commented Mar 31, 2015 at 8:11