---
title: 9-2. Chebyshev approximations - boundary bordering for BVP
tags: Spectral method
---

**Reference**
1. [Spectral method in Matlab](https://people.maths.ox.ac.uk/trefethen/spectral.html) by *Lloyd N. Trefethen*

---

% p13.m - solve linear BVP u_xx = exp(4x), u(-1)=u(1)=0

  N = 16;
  [D,x] = cheb(N);        
  D2 = D^2;              
  D2 = D2(2:N,2:N);                   % boundary conditions
  f = exp(4*x(2:N));           
  u = D2\f;                           % Poisson eq. solved here
  u = [0;u;0];              
  clf, subplot('position',[.1 .4 .8 .5])
  plot(x,u,'.-','markersize',16)
  exact = ( exp(4*x) - sinh(4)*x - cosh(4) )/16; 
  title(['max err = ' num2str(norm(u-exact,inf))],'fontsize',12)

---

% p32.m - solve u_xx = exp(4x), u(-1)=0, u(1)=1 (compare p13.m)

  N = 16;
  [D,x] = cheb(N);
  D2 = D^2;
  D2 = D2(2:N,2:N);
  f = exp(4*x(2:N));           
  u = D2\f;       
  u = [0;u;0] + (x+1)/2;              
  clf
  subplot('position',[.1 .4 .8 .5])
  plot(x,u,'.-','markersize',16)
  xx = -1:.01:1;
  exact = (exp(4*x) - sinh(4)*x - cosh(4))/16 + (x+1)/2;
  title(['max err = ' num2str(norm(u-exact,inf))],'fontsize',12)

---

% p33.m - solve linear BVP u_xx = exp(4x), u'(-1)=u(1)=0

  N = 16;
  [D,x] = cheb(N); D2 = D^2;              
  D2(N+1,:) = D(N+1,:);            % Neumann condition at x = -1
  D2 = D2(2:N+1,2:N+1);                
  f = exp(4*x(2:N));           
  u = D2\[f;0];                       
  u = [0;u];              
  clf, subplot('position',[.1 .4 .8 .5])
  plot(x,u,'.-','markersize',16)
  axis([-1 1 -4 0])
  exact = (exp(4*x) - 4*exp(-4)*(x-1) - exp(4))/16; 
  title(['max err = ' num2str(norm(u-exact,inf))],'fontsize',12)

---

% p34.m - Allen-Cahn eq. u_t = eps*u_xx+u-u^3, u(-1)=-1, u(1)=1
%         (compare p6.m and p32.m)

% Differentiation matrix and initial data:
  N = 20; [D,x] = cheb(N); D2 = D^2;     % use full-size matrix
  D2([1 N+1],:) = zeros(2,N+1);          % for convenience
  eps = 0.01; dt = min([.01,50*N^(-4)/eps]);
  t = 0; v = .53*x + .47*sin(-1.5*pi*x);

% Solve PDE by Euler formula and plot results:
  tmax = 100; tplot = 2; nplots = round(tmax/tplot);
  plotgap = round(tplot/dt); dt = tplot/plotgap;
  plotdata = [v'; zeros(nplots,length(x))]; tdata = t;
  for i = 1:nplots
    for n = 1:plotgap
      t = t+dt; v = v + dt*(eps*D2*(v-x) + v - v.^3);    % Euler
    end
    plotdata(i+1,:) = v; tdata = [tdata; t];
  end
  clf, subplot('position',[.1 .4 .8 .5])
  mesh(x,tdata,plotdata), grid on, axis([-1 1 0 tmax -1 1]),
  view(-60,55), colormap(1e-6*[1 1 1]); xlabel x, ylabel t, zlabel u

---

% p35.m - Allen-Cahn eq. as in p34.m, but with boundary condition
%         imposed explicitly ("method (II)")

% Differentiation matrix and initial data:
  N = 20; [D,x] = cheb(N); D2 = D^2;    
  eps = 0.01; dt = min([.01,50*N^(-4)/eps]);
  t = 0; v = .53*x + .47*sin(-1.5*pi*x);

% Solve PDE by Euler formula and plot results:
  tmax = 100; tplot = 2; nplots = round(tmax/tplot);
  plotgap = round(tplot/dt); dt = tplot/plotgap;
  plotdata = [v'; zeros(nplots,length(x))]; tdata = t;
  for i = 1:nplots
    for n = 1:plotgap
      t = t+dt; v = v + dt*(eps*D2*v + v - v.^3);        % Euler
      v(1) = 1 + sin(t/5)^2;  v(end) = -1;               % BC
    end
    plotdata(i+1,:) = v; tdata = [tdata; t];
  end
  clf, subplot('position',[.1 .4 .8 .5])
  mesh(x,tdata,plotdata), grid on, axis([-1 1 0 tmax -1 2]),
  view(-60,55), colormap(1e-6*[1 1 1]); xlabel x, ylabel t, zlabel u

---

% p38.m - solve u_xxxx = exp(x), u(-1)=u(1)=u'(-1)=u'(1)=0
%         (compare p13.m)

% Construct discrete biharmonic operator:
  N = 15; [D,x] = cheb(N);
  S = diag([0; 1 ./(1-x(2:N).^2); 0]);
  D4 = (diag(1-x.^2)*D^4 - 8*diag(x)*D^3 - 12*D^2)*S;
  D4 = D4(2:N,2:N);

% Solve boundary-value problem and plot result:
  f = exp(x(2:N)); u = D4\f; u = [0;u;0];
  clf, subplot('position',[.1 .4 .8 .5])
  plot(x,u,'.','markersize',16)
  axis([-1 1 -.01 .06]), grid on
  line(x,u)

% Determine exact solution and print maximum error:
  A = [1 -1 1 -1; 0 1 -2 3; 1 1 1 1; 0 1 2 3];
  V = vander(x); V = V(:,end:-1:end-3);
  c = A\exp([-1 -1 1 1]'); exact = exp(x) - V*c;
  title(['max err = ' num2str(norm(u-exact,inf))],'fontsize',12)

---