Plot absolute error vs step-size as a log-log plot (MATLAB)

Question

I have approximated up to $t=2$ with step sizes $h=(1/2)^{n}$ for $n=1, ... , 10$ using a Runge-Kutta (order 4) MATLAB code. Then, I wrote down the absolute error of $y(2)$. Now, I want to plot the log(error) against the log(h), in MATLAB.

I should get a line with slope 4 (I think), because the Runge-Kutta method has order 4 and because $error\approx C(h^{k})$ and so $\log(error)=\log(C)+k\log(h)$, where $k$ is the order of the numerical method. I've tried this several different ways by searching the internet, but nothing seems to be correct. Most recently, I've tried to define a vector $xval$ containing the $h$ step sizes and a vector $yval$ containing the errors. Then I used $loglog(x, y, '-s')$ to obtain a graph, but when I used $polyfit(xval, yval, 1)$, I did not get back a slope of 4.

I'm really confused about this. Any help would be greatly appreciated. Here is my data (step size, error): \begin{align*} &1/2, \qquad 6.99e^{-2}\\ &1/4 \qquad 1.08e^{-2}\\ &1/8 \qquad 1.02e^{-3}\\ &1/16 \qquad 7.2e^{-5}\\ &1/32 \qquad 4.54e^{-6}\\ &1/64 \qquad 2.78e^{-7}\\ &1/128 \qquad 1.71e^{-8}\\ &1/256 \qquad 1.06e^{-9}\\ &1/512 \qquad 6.57e^{-11}\\ &1/1024 \qquad 4.11e^{-12}\\ \end{align*}

Answer 1

I'm using R and obtain

## define new variables:
log_h   = log(h),
log_err = log(err)

## perform a linear fit:
lm(formula = log_err ~ log_h, data = df)

## display the summary of fit:
Residuals:
     Min       1Q   Median       3Q      Max 
-0.76940 -0.13870  0.04283  0.25340  0.37336 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.78235    0.24929   3.138   0.0138 *  
log_h        3.85725    0.05796  66.546 2.89e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3649 on 8 degrees of freedom
Multiple R-squared:  0.9982,    Adjusted R-squared:  0.998 
F-statistic:  4428 on 1 and 8 DF,  p-value: 2.893e-12

So the estimated slope is $3.85 \pm 0.05$, which is pretty close to 4. The plot shows a straight line. The residuals are OK -- not shown here.

In Matlab one would do something like

p = polyfit(log_h, log_err, 1)