% % chap2fig2.m % % second order finite difference method for the bvp % u'' = f, u'(ax)=sigma, u(bx)=beta % Using 5-pt differences on an arbitrary grid. % % Produces Figure 2.2 % From http://www.amath.washington.edu/~rjl/fdmbook/chapter2 (2007) ax = 0; bx = 1; sigma = 0; beta = 3; f = @(x) exp(x); % right hand side function utrue = @(x) exp(x) + (sigma-exp(ax))*(x - bx) + beta - exp(bx); % true soln % Solve the problem for ntest different grid sizes to test convergence: m1vals = [10 20 40 80]; ntest = length(m1vals); hvals = zeros(ntest,1); % to hold h values E = zeros(ntest,1); % to hold errors u11 = zeros(11,2); % to store 11-point solution for plotting for approach=1:2 for jtest=1:ntest m1 = m1vals(jtest); m2 = m1 + 1; m = m1 - 1; h = (bx-ax)/m1; hvals(jtest) = h; x = linspace(ax,bx,m2)'; % m+2 points counting the boundaries % set up matrix A (using sparse matrix storage): e = ones(m2,1); A = 1/h^2 * spdiags([e -2*e e], [-1 0 1], m2, m2); % fix up values in first and last row for BCs: A(1,1) = -1/h; A(1,2) = 1/h; A(m2,m1) = 0; A(m2,m2) = 1; % Right hand side: F = f(x); F(1) = sigma + h/2 * F(1) * (approach-1); % Approach 1 or 2 from Section 2.12 F(m2) = beta; % solve linear system: U = A\F; % compute errors: uhat = utrue(x); err = U - uhat; E(jtest,approach) = max(abs(err)); if jtest==1 % save for plotting: x1 = x; u11(:,approach) = U; end end end % plot solutions: figure(1) clf xfine = linspace(ax,bx,101); ufine = utrue(xfine); plot(xfine,ufine) hold on plot(x1,u11(:,1),'+') plot(x1,u11(:,2),'o') hold off axis([-.1 1.1 2.2 3.2]) set(gca,'fontsize',15) % log-log plot of errors: figure(2) clf loglog(hvals,E(:,1),'+-') hold on loglog(hvals,E(:,2),'o-') hold off axis([5e-3 .2 1e-5 1]) set(gca,'fontsize',15)