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

Question

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$:

Answer 1

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$.

Answer 2

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.

Answer 3

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$$