Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

15 Jul 2011 (Updated )

The MATLAB package Chebpack solves specific problems for differential or integral equations.

Boussinesq.m
```function [lam,phi,phip,phipp,t,C,mass]=Boussinesq(n,numeigval,dom,dt,tfin)
%
% From: A. Mohebbi,Z. Asgari, Efficient numerical algorithms for the
% solution of good Boussinesq equation in water wave propagation,
% Computer Physics Communications 182 (2011) 24642470
%
% u_tt+u_xxxx-u_xx-(u^2)_xx=0, [-50,50]
% u(x,0)=uex(x,0), u_t(x,0)=uex_t(x,0)
% u(-50,t)=u(50,t)=u_x(-50,t)=u_x(50,t)=0
%
% Exact solution: uex(x,t)=-{Asech^2(sqrt(A/6)(x-ct-x0))+b+1/2}
%                 A=0.369, b=-0.5, x0=0, c=sqrt(-2(b+A/3))
% Call: [lam,phi,phip,phipp,t,C,mass]=Boussinesq(256,100,[-50,50],0.1,20);
% Boussinesq problem
%
A=0.369;b=-0.5;x0=0;c=sqrt(-2*(b+A/3));
phi=zeros(n,numeigval);phip=phi;phipp=phi;
disp('Linear part');tic;
E=speye(n);D=deriv(n,dom);X=mult(n,dom);
[xleg,wleg]=pd(n,dom,3);Tleg=cpv(n,xleg,dom);rleg=ones(n,1);
F=(dom(1)*E-X)^2*(dom(2)*E-X)^2;
AM=F*(D^4-D^2)*F;B=F*F;AA=AM(1:n-4,1:n-4);BB=B(1:n-4,1:n-4);
[V,L]=eig(full(AA),full(BB));[LL,ind]=sort(diag(L),'ascend');VV=V(:,ind);
lam=LL(1:numeigval);
for s=1:numeigval
phi(:,s)=Tleg*F*[VV(:,s);0;0;0;0];nphi=sqrt(wleg*phi(:,s).^2);
phi(:,s)=phi(:,s)/nphi;
phip(:,s)=Tleg*D*F*[VV(:,s);0;0;0;0]/nphi;
phipp(:,s)=Tleg*D^2*F*[VV(:,s);0;0;0;0]/nphi;
end
toc;
disp(['Calculated ',num2str(numeigval),' eigenvalues and eigenfunctions']);
I=wip(phi,phi,rleg,wleg);orthogonal=max(max(abs(I-eye(numeigval))));
disp(['test orthogonality: ',num2str(orthogonal)]);tic;
u0=ICB(xleg);c0=[wip(phi,u0(:,1),rleg,wleg);wip(phi,u0(:,2),rleg,wleg)];
opts=odeset('RelTol',1.e-9,'AbsTol',1.e-9);
[t,C]=ode15s(@Bfun,0:dt:tfin,c0,opts);lt=length(t);
solv=[phi*C(:,1:numeigval)';phip*C(:,numeigval+1:2*numeigval)'];toc;
err=zeros(lt,1);mass=zeros(lt,1);
for k=1:lt
uex=-(A*sech(sqrt(A/6)*(xleg-c*t(k)-x0)).^2+b+0.5);
err(k)=norm(uex-solv(1:n,k));
mass(k)=wleg*solv(1:n,k);
end
figure(1)
surf((xleg),t,-(solv(1:n,:))'), shading interp, lighting phong, axis tight
view([30 70]), colormap(autumn); set(gca,'zlim',[0 0.5])
light('color',[1 1 0],'position',[-1,2,2])
material([0.30 0.60 0.60 40.00 1.00]);
xlabel('x');ylabel('t');zlabel('-y(x,t)');
figure(2);semilogy(t,err);title('error');grid;xlabel('t');ylabel('err');
figure(3);plot(t,mass-mass(1));title('invariant');
grid;xlabel('t');ylabel('mass-mass(0)');
function U=ICB(x)
A=0.369;b=-0.5;x0=0;c=sqrt(-2*(b+A/3));
U(:,1)=-(A*sech(sqrt(A/6)*(x-x0)).^2+b+0.5);% function
U(:,2)= 2*A*sech(sqrt(A/6)*(x-x0)).^2.*tanh(sqrt(A/6)*(x-x0))*...
(-c*sqrt(A/6));% the derivative
end
function dcdt=Bfun(~,cc)
u=phi*cc(1:numeigval);up=phip*cc(1:numeigval);upp=phipp*cc(1:numeigval);
Nu=2*up.^2+2*u.*upp;
coef=wip(phi,Nu,rleg,wleg);
dcdt=[cc(numeigval+1:2*numeigval);-diag(lam)*cc(1:numeigval)+coef];
end
end```