Computing and determining the convergence rate

Question

Write a MATLAB code to compute and determine the convergence rate of :

(exp(h)-(1+h+1/2*h^2))/h ‍‍‍‍‍‍ ‍‍‍‍‍‍ with h=1/2, 1/2^2,..., 1/2^10

My code was:

h0=(0.5)^i;
TOL=10^(-8);
N=10;
i=1;
flag=0;
table=zeros(30,1);
table(1)=h0

while i < N
    h=(exp(h0)-(1+h0+0.5*h0^2))/h0;
    table (i+1)=h;
    if abs(h-h0)< TOL
        flag=1;
        break;
    end
    i=i+1;
    h0=h;
end

if flag==1
    h
else
    error('failed');
end

The answer I received does not make quite any sense. Please help.

One bug I see in your code is that you're using `h0=(0.5)^i` which means `h0=(0.5)^√-1;` which means `h0=9.7656e-04`. I think you wanted to use `h0=0.5`. And Questions seeking debugging help (*"why isn't this code working?"*) must include the desired behavior. Tell us what your expected answer is! — Sardar Usama, Sep 12 '16 at 06:01

NLindros · Answer 1 · 2016-09-12T06:38:25.830

The main problem is in your tolerance check

if abs(h-h0) < TOL

If your expression approaches its limit fast enough, h can become 0 as h0 is larger than the tolerance. If so, the criteria isn't met and the loop continues. Next iteration h0 is 0 and h will be evaluated as NaN (since divition with 0 is bad)

If you, like in this case, expect a convergence towards 0, you could instead check

if h > TOL

or you could also add a NaN-check

 if abs(h-h0) < TOL || isnan(h)

Apart from that, there are a few problems with your code.

First you initiate h0 using i (the imaginary number), you probably intented to use i = 1, but that line is below in you code.

You use a whileloop but in your case, as you intend to increment i, a forloop would be just as good.

The use of both the variable table, h and h0 is unnecessary. Make a single result vector, initialized with your h0 at the first index - see example below.

TOL = 1e-8; % Tolerance
N = 10; % Max number of iterations

% Value vector
h = zeros(N+1,1);

% Init value
h(1) = (0.5)^1;

for k = 1:N
    h(k+1) = (exp(h(k)) - (1 + h(k) + 0.5*h(k)^2))/h(k);
    if isnan(h(k+1)) || abs(diff(h(k + [0 1]))) < TOL
        N = k;
        break
    end
end
% Crop vector
h = h(1:N);

% Display result
fprintf('Converged after %d iterations\n', N)

% Plot (with logarithmic y-xis)
semilogy(1:N, h,'*-')

Computing and determining the convergence rate

1 Answers1