1
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ does the solution $y$ satisfy some ODE, possibly of the form $y'=f(t,y)$? $\endgroup$ – Ellya Mar 30 '15 at 23:53
  • $\begingroup$ @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$? $\endgroup$ – evinda Mar 31 '15 at 11:31
  • $\begingroup$ 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. $\endgroup$ – LutzL Mar 31 '15 at 11:44
  • $\begingroup$ 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. $\endgroup$ – Ross Millikan Mar 31 '15 at 14:12
  • 1
    $\begingroup$ 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.) $\endgroup$ – LutzL Apr 1 '15 at 7:52
2
+50
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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 )}$$ ? $\endgroup$ – evinda Apr 1 '15 at 22:14
  • $\begingroup$ 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|$. $\endgroup$ – LutzL Apr 2 '15 at 7:30
  • $\begingroup$ I edited my post... Could you take a look at it? $\endgroup$ – evinda Apr 2 '15 at 8:32
  • $\begingroup$ @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. $\endgroup$ – LutzL Apr 2 '15 at 9:02
  • 1
    $\begingroup$ 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. $\endgroup$ – LutzL Apr 2 '15 at 18:40
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ To do numeric integration with $N_1,N_2$ steps do we not have to take a specific $\rho$ ? $\endgroup$ – evinda Mar 31 '15 at 11:23
  • $\begingroup$ 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. $\endgroup$ – Ross Millikan Mar 31 '15 at 14:15
  • $\begingroup$ math.stackexchange.com/questions/1197371/… Could you take a look at my question? $\endgroup$ – evinda Apr 30 '15 at 17:08
0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ @LutzL cheers edited :) $\endgroup$ – Ellya Mar 31 '15 at 7:27
  • $\begingroup$ @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 $\endgroup$ – Ellya Mar 31 '15 at 7:30
  • $\begingroup$ 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. $\endgroup$ – LutzL Mar 31 '15 at 7:36
  • $\begingroup$ 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$. $\endgroup$ – LutzL Mar 31 '15 at 8:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.