% SCALAR_FIG1 Figures 1-3 from scalar.tex function scalar_fig1 L = 30; % length of each of the halves %xik = 0; xik = -1; % where to evaluate local error xistar = 0; % where to evaluate wedge product xim = xistar - L; xip = xistar + L; hs = [ 0.2 0.1 ]; hsmall = 0.02; las = logspace(0,8,200)*1i; c = -5/6*sqrt(6); kas = sqrt(c^2 + 4*(las-1)); %hs = 0.1; las = 1i; %hsmall = 1; %las = 1000i; ol = zeros(length(hs), length(las), 2); for m = 1:length(las) la = las(m); mu1 = 1/2*(-c+sqrt(c^2+4*(la+1))); mu2 = 1/2*(-c-sqrt(c^2+4*(la+1))); B = [ 1, 1; mu1, mu2 ]; Bi = inv(B); y0 = B * [1;0]; yk = solve_gl4(xim, y0, xik, hsmall, la, mu1); yzerom = solve_gl4(xik, yk, xistar, hsmall, la, mu1); mu1p = 1/2*(-c+sqrt(c^2+4*(la-1))); mu2p = 1/2*(-c-sqrt(c^2+4*(la-1))); Bp = [ 1, 1; mu1p, mu2p ]; y0p = Bp * [0;1]; yzerop = solve_gl4(xip, y0p, xistar, -hsmall, la, mu2p); Dex = yzerom(1)*yzerop(2) - yzerom(2)*yzerop(1); %foobar; for n = 1:length(hs) % Magnus-4 h = hs(n); y1 = solve_gl4(xik, yk, xik+h, hsmall, la, mu1); y2 = solve_m4(xik, yk, xik+h, h, la, mu1); ol(n,m,:) = Bi*(y2-y1); y3 = solve_m4(xim, y0, xik, h, la, mu1); og(n,m,:) = Bi*(y3-yk); yzerom = solve_m4(xik, y3, xistar, h, la, mu1); yzerop = solve_m4(xip, y0p, xistar, -h, la, mu2p); ed(n,m) = yzerom(1)*yzerop(2) - yzerom(2)*yzerop(1) - Dex; % Gauss-Legendre y2 = solve_gl4(xik, yk, xik+h, h, la, mu1); olgl(n,m,:) = Bi*(y2-y1); y3 = solve_gl4(xim, y0, xik, h, la, mu1); oggl(n,m,:) = Bi*(y3-yk); yzerom = solve_gl4(xim, y0, xistar, h, la, mu1); yzerop = solve_gl4(xip, y0p, xistar, -h, la, mu2p); edgl(n,m) = yzerom(1)*yzerop(2) - yzerom(2)*yzerop(1) - Dex; % Exponential midpoint yzerom = solve_m2(xim, y0, xistar, h, la, mu1); yzerop = solve_m2(xip, y0p, xistar, -h, la, mu2p); edm2(n,m) = yzerom(1)*yzerop(2) - yzerom(2)*yzerop(1) - Dex; end end; for n = 1:length(hs) % formulas for xik = -1 ole(n,:,1) = 0.00038544 * hs(n)^5 ./ kas; ole(n,:,2) = 0.019607 * hs(n)^2 ./ kas; oge(n,:,1) = 0.0014950 * hs(n)^4 ./ kas; oge(n,:,2) = 0.019607 * hs(n)^2 ./ kas; ede(n,:) = -0.002268 * hs(n)^4; olegl(n,:,1) = -0.0000010183 * hs(n)^5 ./ kas; olegl(n,:,2) = 0.00074113 * hs(n)^3 ./ (kas.^2); ogegl(n,:,1) = -0.0000061761 * hs(n)^4 ./ kas; ogegl(n,:,2) = 0.00074113 * hs(n)^3 ./ (kas.^2); ogegl2(n,:,1) = -0.000046537 * hs(n)^4 ./ (kas.^2); edgle(n,:) = log(2)/51840 * hs(n)^4 ./ sqrt(las); edm2e(n,:) = log(2)/24 * hs(n)^2 ./ sqrt(las); end; figure(1); clf; subplot(2,2,1); loglog(imag(las), abs(ol(1,:,1)), '-', ... imag(las), abs(ole(1,:,1)), '--', ... imag(las), abs(ol(2,:,1)), '-', ... imag(las), abs(ole(2,:,1)), '--'); axis([1e0 1e8 1e-13 1e-7]); text(4e+5, 2e-10, '2'); text(4e+5, 7e-12, '1'); title('First component of local error'); xlabel('|\lambda|'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|[{\itL_k}]_1|'); set(gca, 'YTick', [1e-13 1e-10 1e-7]); subplot(2,2,2); loglog(imag(las), abs(ol(1,:,2)), '-', ... imag(las), abs(ole(1,:,2)), '--', ... imag(las), abs(ol(2,:,2)), '-', ... imag(las), abs(ole(2,:,2)), '--'); axis([1e0 1e8 1e-9 1e-3]); text(1e6, 9e-7, '2'); text(1e6, 3e-8, '1'); title('Second component of local error'); xlabel('|\lambda|'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|[{\itL_k}]_2|'); set(gca, 'YTick', [1e-9 1e-6 1e-3]); subplot(2,2,3); loglog(imag(las), abs(og(1,:,1)), '-', ... imag(las), abs(oge(1,:,1)), '--', ... imag(las), abs(og(2,:,1)), '-', ... imag(las), abs(oge(2,:,1)), '--'); axis([1e0 1e8 1e-11 1e-5]); text(3e+4, 1.5e-8, '2'); text(3e+4, 8e-10, '1'); title('First component of global error'); xlabel('|\lambda|'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|[{\itE_k}]_1|'); set(gca, 'YTick', [1e-11 1e-8 1e-5]); subplot(2,2,4); loglog(imag(las), abs(og(1,:,2)), '-', ... imag(las), abs(oge(1,:,2)), '--', ... imag(las), abs(og(2,:,2)), '-', ... imag(las), abs(oge(2,:,2)), '--'); axis([1e0 1e8 1e-9 1e-3]); text(1e6, 9e-7, '2'); text(1e6, 3e-8, '1'); title('Second component of global error'); xlabel('|\lambda|'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|[{\itE_k}]_2|'); set(gca, 'YTick', [1e-9 1e-6 1e-3]); print -deps scalar_fig1.eps; figure(2); clf; subplot(2,2,1); loglog(imag(las), abs(ed(1,:)), '-', ... imag(las([1 end])), abs(ede(1,[1 end])), '--'); axis([1e0 1e8 1e-8 1e-4]); title('Error in the Evans function for {\ith} = 0.2'); xlabel('\lambda'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|{\itE_D}|'); subplot(2,2,2); loglog(imag(las), abs(ed(2,:)), '-', ... imag(las([1 end])), abs(ede(2,[1 end])), '--'); axis([1e0 1e8 1e-8 1e-4]); title('Error in the Evans function for {\ith} = 0.1'); xlabel('\lambda'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|{\itE_D}|'); print -deps scalar_fig2.eps; figure(3); clf; subplot(2,2,1); loglog(imag(las), abs(edm2(1,:)), '-', ... imag(las), abs(edm2(2,:)), '-');%, ... %imag(las([1 end])), abs(edm2e(1,[1 end])), '--', ... %imag(las([1 end])), abs(edm2e(2,[1 end])), '--'); axis([1e0 1e8 1e-9 1e-4]); text(4e4, 6e-9, '2'); text(6e4, 2.5e-8, '1'); title('Error in {\itD}(\lambda) for exponential midpoint'); xlabel('\lambda'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|{\itE_D}|'); subplot(2,2,2); loglog(imag(las), abs(edgl(1,:)), '-', ... imag(las), abs(edgl(2,:)), '-', ... imag(las), 1e-3 * 0.1^4 ./ abs(las), ':', ... imag(las), 1e-3 * 0.2^4 ./ abs(las), ':'); axis([1e0 1e8 3e-13 3e-7]); text(1e2, 8e-11, '2'); text(4, 1e-7, '1'); title('Error in {\itD}(\lambda) for Gauss-Legendre'); xlabel('\lambda'); set(gca, 'XTick', [1 1e4 1e8]); ylabel('|{\itE_D}|'); print -deps scalar_fig3.eps; % figure(4); clf; % subplot(2,2,1); % loglog(imag(las), abs(olgl(1,:,1)), '-', ... % imag(las), abs(olegl(1,:,1)), '--', ... % imag(las), abs(olgl(2,:,1)), '-', ... % imag(las), abs(olegl(2,:,1)), '--'); % axis([1e0 1e8 1e-13 1e-7]); % text(4e+5, 2e-10, '1'); % text(4e+5, 7e-12, '2'); % title('First component of local error'); % xlabel('|\lambda|'); % set(gca, 'XTick', [1 1e4 1e8]); % ylabel('|[{\itL_k}]_1|'); % set(gca, 'YTick', [1e-13 1e-10 1e-7]); % subplot(2,2,2); % loglog(imag(las), abs(olgl(1,:,2)), '-', ... % imag(las), abs(olegl(1,:,2)), '--', ... % imag(las), abs(olgl(2,:,2)), '-', ... % imag(las), abs(olegl(2,:,2)), '--'); % axis([1e0 1e8 1e-15 1e-8]); % text(1e6, 9e-7, '1'); % text(1e6, 3e-8, '2'); % title('Second component of local error'); % xlabel('|\lambda|'); % set(gca, 'XTick', [1 1e4 1e8]); % ylabel('|[{\itL_k}]_2|'); % set(gca, 'YTick', [1e-15 1e-10 1e-8]); % % subplot(2,2,1); % loglog(imag(las), abs(oggl(1,:,1)), '-', ... % imag(las), abs(ogegl(1,:,1)), '--', ... % imag(las), abs(oggl(2,:,1)), '-', ... % imag(las), abs(ogegl(2,:,1)), '--'); % axis([1e0 1e8 1e-13 1e-7]); % % text(3e+4, 1.5e-8, '1'); % % text(3e+4, 8e-10, '2'); % title('First component of global error'); % xlabel('|\lambda|'); % set(gca, 'XTick', [1 1e4 1e8]); % ylabel('|[{\itE_k}]_1|'); % set(gca, 'YTick', [1e-13 1e-10 1e-7]); % % subplot(2,2,4); % loglog(imag(las), abs(oggl(1,:,2)), '-', ... % imag(las), abs(ogegl(1,:,2)), '--', ... % imag(las), abs(oggl(2,:,2)), '-', ... % imag(las), abs(ogegl(2,:,2)), '--'); % axis([1e0 1e8 1e-13 1e-7]); % % text(1e6, 9e-7, '1'); % % text(1e6, 3e-8, '2'); % title('Second component of global error'); % xlabel('|\lambda|'); % set(gca, 'XTick', [1 1e4 1e8]); % ylabel('|[{\itE_k}]_2|'); % set(gca, 'YTick', [1e-13 1e-10 1e-7]); % % subplot(2,2,3); % loglog(imag(las), abs(oggl(1,:,1)-ogegl(1,:,1)), '-', ... % imag(las), abs(ogegl2(1,:,1)), '--', ... % imag(las), abs(oggl(2,:,1)-ogegl(2,:,1)), '-', ... % imag(las), abs(ogegl2(2,:,1)), '--'); % axis([1e0 1e8 1e-15 1e-7]); % % text(1e6, 9e-7, '1'); % % text(1e6, 3e-8, '2'); % title('Global error minus estimate (1st comp.)'); % xlabel('|\lambda|'); % set(gca, 'XTick', [1 1e4 1e8]); % ylabel('|[{\itE_k}]_1 - est_1|'); % set(gca, 'YTick', [1e-15 1e-11 1e-7]); % figure(4); clf; % loglog(imag(las), abs(ed(1,:)), '-', ... % imag(las), abs(ed(2,:)), '-'); % title('Error in {\itD}(\lambda) for Magnus-4'); % xlabel('\lambda'); % ylabel('|{\itE_D}|'); % print -deps scalar.fig5.eps; % % figure(4); clf; % loglog(imag(las), abs(edm2(1,:)), '-', ... % imag(las), abs(edm2(2,:)), '-'); % title('Error in {\itD}(\lambda) for exponential midpoint'); % xlabel('\lambda'); % ylabel('|{\itE_D}|'); % print -deps scalar.fig6.eps; % % figure(4); clf; % loglog(imag(las), abs(edgl(1,:)), '-', ... % imag(las), abs(edgl(2,:)), '-'); % title('Error in {\itD}(\lambda) for Gauss-Legendre'); % xlabel('\lambda'); % ylabel('|{\itE_D}|'); % print -deps scalar.fig7.eps; end function yk = solve_gl4(xi0, y0, xik, h_param, la, mu) nofsteps = round((xik-xi0) / h_param); if nofsteps ~= 0 h = (xik-xi0) / nofsteps; else h = 0; end xis = linspace(xi0, xik, nofsteps+1); A1 = fisher_A(xis + (1/2-1/6*sqrt(3))*h, la, mu); A2 = fisher_A(xis + (1/2+1/6*sqrt(3))*h, la, mu); yk = y0; %ys = [y0]; for k = 1:nofsteps A1k = A1(:,:,k); A2k = A2(:,:,k); A1k*yk ./ yk; U = [ eye(2)-h/4*A1k, -(1/4-sqrt(3)/6)*h*A1k -(1/4+sqrt(3)/6)*h*A2k, eye(2)-h/4*A2k ]; kk = U \ [ A1k*yk; A2k*yk ]; yk = yk + h/2 * (kk(1:2) + kk(3:4)); %ys = [ys yk]; end %plot(xis, abs(ys), '-o'); end function yk = solve_m4(xi0, y0, xik, h_param, la, mu) nofsteps = round((xik-xi0) / h_param); if nofsteps ~= 0 h = (xik-xi0) / nofsteps; else h = 0; end xis = linspace(xi0, xik, nofsteps+1); A1 = fisher_A(xis + (1/2-1/6*sqrt(3))*h, la, mu); A2 = fisher_A(xis + (1/2+1/6*sqrt(3))*h, la, mu); yk = y0; for k = 1:nofsteps A1k = A1(:,:,k); A2k = A2(:,:,k); om = h/2 * (A1k+A2k) - sqrt(3)/12*h^2 * (A1k*A2k-A2k*A1k); % yk = expm(om) * yk; yk = expmdemo3(om) * yk; % yk = expmanal(om) * yk; end end function yk = solve_m2(xi0, y0, xik, h_param, la, mu) nofsteps = round((xik-xi0) / h_param); if nofsteps ~= 0 h = (xik-xi0) / nofsteps; else h = 0; end xis = linspace(xi0, xik, nofsteps+1); A1 = fisher_A(xis + 1/2*h, la, mu); yk = y0; for k = 1:nofsteps A1k = A1(:,:,k); yk = expmdemo3(h*A1k) * yk; end end function res = fisher_A(xis, la, mu); res = zeros(2, 2, length(xis)); res(1,1,:) = -mu; res(1,2,:) = 1; res(2,1,:) = la - 1 + 2 ./ (1+exp(xis/sqrt(6))).^2; res(2,2,:) = 5/6*sqrt(6)-mu; end function res = expmanal(A) % Analytic formula for exponential of 2x2 matrices, from % http://mathworld.wolfram.com/MatrixExponential.html % Does not work if de = 0. a = A(1,1); b = A(1,2); c = A(2,1); d = A(2,2); de = sqrt((a-d)^2 + 4*b*c); res = exp((a+d)/2) / de * ... [ de*cosh(de/2)+(a-d)*sinh(de/2) 2*b*sinh(de/2) 2*c*sinh(de/2) de*cosh(de/2)+(d-a)*sinh(de/2) ]; end