1function [err,rate] = selfConvergence(ug,h,M)
3if nargin == 0, test_selfConvergence; return; end
5if(M<3)
6 error('need at least 3 resolutions to perform self-convergence study');
7end
9m = M-2;
11Dmp1 = ug{m+1} - ug{m};
12Dmp2 = ug{m+2} - ug{m+1};
14NDmp1 = maxNorm(Dmp1);
15NDmp2 = maxNorm(Dmp2);
16ratio = NDmp1/NDmp2;
17r = h(m)/h(m+1);
18rate = log(ratio)/log(r);
20C = Dmp2/abs((h(m+1)^(rate)) - (h(m+2)^(rate)));
22err = zeros(1,M);
23for k=1:M
24 uerr = C*(h(k)^(rate));
25 err(k) = maxNorm(uerr);
26end
28fprintf('self convergence rate = %1.2e ', rate);
29fmt=['errors =' repmat(' %1.2e;',1,numel(err)) '\n'];
30fprintf(fmt,err);
31end
33%%% test function
34function test_selfConvergence
36ue = @(x,t) cos(x + 3).*exp(x + 1).*log(1 + x).*sin(2*t - 2).*exp(-t);
38h0 = 1e-3;
39h(1) = h0;
40h(2) = h0/2;
41h(3) = h0/4;
42h(4) = h0/8;
44x = -pi:h0:pi;
45t = linspace(0,1,length(x));
47rate = 1.5;
48numRes = 4;
49ug = cell(numRes,1);
50C = (1 + (100-1).*rand(size(x)));
51for k = 1:numRes
52 ug{k} = ue(x,t) + C*(h(k)^(rate)) + C*(h(k)^(rate+3));
53end
55selfConvergence(ug,h,numRes);
57end