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