% % bvpspectral.m % pseudospectral finite difference method for the bvp % u''(x) = f(x), u'(ax)=sigma, u(bx)=beta % The use of Chebyshev grid points gives a very accurate method. % With a uniform grid it is quite accurate for small m, but unstable as % m increases. % % Different BCs can be specified by changing the first and/or last rows of % A and F. % % From http://www.amath.washington.edu/~rjl/fdmbook/chapter2 (2007) ax = 0; bx = 3; sigma = -5; % Neumann boundary condition at ax beta = 3; % Dirichlet boundary condtion at bx f = @(x) exp(x); % right hand side function utrue = @(x) exp(x) + (sigma-exp(ax))*(x - bx) + beta - exp(bx); % true soln % true solution on fine grid for plotting: xfine = linspace(ax,bx,101); ufine = utrue(xfine); % Solve the problem for ntest different grid sizes to test convergence: m1vals = [5 10 15 20]; % note: fewer points than used for other methods ntest = length(m1vals); hvals = zeros(ntest,1); % to hold h values E = zeros(ntest,1); % to hold errors for jtest=1:ntest m1 = m1vals(jtest); m2 = m1 + 1; m = m1 - 1; % number of interior grid points hvals(jtest) = (bx-ax)/m1; % average grid spacing, for convergence tests % set grid points: gridchoice = 'chebyshev'; x = xgrid(ax,bx,m,gridchoice); % set up matrix A: A = zeros(m2,m2); % initialize (A is dense) % first row for Neumann BC at ax: A(1,:) = fdcoeffF(1, x(1), x); % interior rows: for i=2:m1 A(i,:) = fdcoeffF(2, x(i), x); end % last row for Dirichlet BC at bx: A(m2,:) = fdcoeffF(0,x(m2),x); % Right hand side: F = f(x); F(1) = sigma; F(m2) = beta; % solve linear system: U = A\F; % compute error at grid points: uhat = utrue(x); err = U - uhat; E(jtest) = max(abs(err)); disp(' ') disp(sprintf('Error with %i points is %9.5e',m2,E(jtest))) clf plot(x,U,'o') % plot computed solution title(sprintf('Computed solution with %i grid points',m2)); hold on plot(xfine,ufine) % plot true solution hold off % pause to see this plot: drawnow input('Hit for next plot '); end error_table(hvals, E); % print tables of errors and ratios error_loglog(hvals, E); % produce log-log plot of errors and least squares fit