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.

[xm,yvm]=test3(n,dom,kind,p,ep,xim,yim)
```function [xm,yvm]=test3(n,dom,kind,p,ep,xim,yim)
% p: order of the equation, ep: the continuation parameter
% Example:  y^(4) = ep(y'y''-yy^(3)), x in [0,1], ep=10000
%           y(0)=y'(0)=0, y(1)=1, y'(1)=0
% From: J. Cash, Test problems for BVP_software, test problem 32,
%       http://www2.imperial.ac.uk/~jcash/BVP_software/PROBLEMS.PDF
% call:  [x,yv]=test3(32,[0 1],2,4,1);
% then   [x,yv]=test3(32,[0 0.2 0.4 1],2,4,10,x,yv);
%        [x,yv]=test3(32,[0 0.2 0.4 1],2,4,100,x,yv);
%        [x,yv]=test3(32,[0 0.2 0.4 1],2,4,1000,x,yv);
%        [x,yv]=test3(32,[0 0.2 0.4 1],2,4,10000,x,yv);
%
xs=pd(n,[-1,1],kind);Xs=mult(n,[-1,1]);Ds=deriv(n,[-1,1]);
m=length(dom);T=cpv(n,[-1,1],[-1,1]);Tl=T(1,:);Tr=T(2,:);
xm=xs*(dom(2:m)-dom(1:m-1))/2+ones(n,1)*(dom(2:m)+dom(1:m-1))/2;
x=reshape(xm,n*(m-1),1);arg=nargin;
AA=spalloc(n*(m-1),n*(m-1),100*n*(m-1));
INIT;% first iteration
dy=1;cont=1;
while dy(cont) > 1.e-12
if cont>100, break, end;
cont=cont+1;
[x,~,sol]=linear_split(n,dom,kind,p,ep,yv,y);
dy(cont)=norm(sol-y);y=sol;
ym=reshape(y,n,m-1);yvm=t2x(ym,kind);
yv=reshape(yvm,n*(m-1),1);
end;
xm=reshape(x,n,m-1);
figure(1);% out
subplot(2,1,1);semilogy(xm,abs(ym),'.');grid;
title('The absolute value of the coefficients of the solution');
subplot(2,1,2);plot(xm,yvm,dom,zeros(m,1),'.r');grid;
xlabel('x');ylabel('y(x)');title('The numerical solution');
figure(2);
plot(dy,'.');
title('The Newton corrections');ylabel('dy');grid;
myOUT;

function INIT
% starting iterations
%
yv=zeros(n*(m-1),1);I=zeros(n*(m-1),1);
if arg<6 % first iteration
yvm=myINIT(xm);ym=x2t(yvm,kind);
yv=reshape(yvm,n*(m-1),1);y=reshape(ym,n*(m-1),1);
else [ni,mi]=size(xim);
for i=1:mi
I(x>=xim(1,i)&x<=xim(ni,i))=i;
yv(I==i)=barycheb(x(I==i),yim(:,i),xim(:,i),kind);
end
yvm=reshape(yv,n,m-1);ym=x2t(yvm,kind);
y=reshape(ym,n*(m-1),1);
end
end

function [x,solv,sol]=linear_split(n,dom,kind,p,ep,yv,y)
% solves the linearized equation
% yv: the old solution in physical form
% y : the old solution in spectral form
%
AA=spalloc(n*(m-1),n*(m-1),100*n*(m-1));
bb=zeros(n*(m-1),1);xx=[];
for k=1:m-1 % the large matrix of the problem
pdom=[dom(k),dom(k+1)]; % for each subinterval
l=(pdom(2)-pdom(1))/2;med=(pdom(2)+pdom(1))/2;
x=l*xs+med;D=Ds/l;X=l*Xs+med*speye(n);xx=[xx,x];
yvk=yv((k-1)*n+1:k*n);yk=y((k-1)*n+1:k*n);
dr=df(x,yvk,yk);% calculates df/dy, df/dy',...
[nl,ml]=size(dr);DF=zeros(n);
for j=1:ml
DF=DF+fact(x2t(dr(:,j),kind),n)*D^(j-1);
end
myDE;A=AL-DF;r=f(x,yvk,yk);b=-DF*yk+x2t(r,kind);% r.h.s.
E=zeros(m-1);E(k,k)=1;A(n-p+1:n,:)=0;AA=AA+kron(E,A);
EE=zeros(m-1,1);EE(k)=1;b(n-p+1:n)=0;bb=bb+kron(EE,b);
end;
myBC;
for j=2:m-1 % the smoothness conditions
fl=(dom(j)-dom(j-1))/2;fr=(dom(j+1)-dom(j))/2;
jnp=j*n-p;j2n=(j-2)*n+1:(j-1)*n;j1n=(j-1)*n+1:j*n;
AA(jnp+1,j2n)=Tr;
AA(jnp+1,j1n)=-Tl;bb(jnp+1)=0;
for k=2:p
AA(jnp+k,j2n)=AA(jnp+k-1,j2n)*Ds/fl;
AA(jnp+k,j1n)=AA(jnp+k-1,j1n)*Ds/fr;
bb(jnp+k)=0;
end
end
sol=AA\bb;solnum=[];
for i=1:m-1
solnum = [solnum,t2x(sol((i-1)*n+1:i*n),kind)];
end
x=reshape(xx,(m-1)*n,1);solv=reshape(solnum,(m-1)*n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          The following functions must be written by the user            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
% linear part of the equation Ly=f(x,y) in spectral form
%
AL=D^4;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
% y(0)=..;y(1)=...
%
AA(n-3,1:n)=Tl;bb(n-3)=0;AA(n-2,1:n)=Tl*Ds*2/(dom(2)-dom(1));bb(n-2)=0;
AA(n-1,(m-2)*n+1:(m-1)*n)=Tr;bb(n-1)=1;
AA(n,(m-2)*n+1:(m-1)*n)=Tr*Ds*2/(dom(m)-dom(m-1));bb(n)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rez=f(x,yv,y)
% nonlinear part of the equation Ly=f(x,y) in physical form
% x: the grid,
% yv: the solution in physical form,
% y: the solution in spectral form
% rez: the values of f
%
yp=t2x(D*y,kind);ys=t2x(D*(D*y),kind);yt=t2x(D*(D*(D*y)),kind);
rez=ep*(yp.*ys-yv.*yt);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drez=df(x,yv,y)
% calculates f_y,f_y',...in drez
%
yp=t2x(D*y,kind);ys=t2x(D*(D*y),kind);yt=t2x(D*(D*(D*y)),kind);
drez(:,1)=-ep*yt;drez(:,2)=ep*ys;drez(:,3)=ep*yp;drez(:,4)=-ep*yv;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yv=myINIT(x)
% starting approximation in the physical space
yv=3*x.^2/4-x.^3/2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
figure(3);
plot(xm,t2x((Ds*ym)*diag(2./(dom(2:m)-dom(1:m-1))),kind),dom,zeros(m,1),'.r');
grid;xlabel('x');ylabel('y''(x)');
title('The derivative of the numerical solution');
%    figure(4);
%    spy(AA);title('spy(AA)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end```