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*}
$\endgroup$ 41 Answer
$\begingroup$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-12So 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)
