Code covered by the BSD License  

Highlights from
Chebpack

image thumbnail

Chebpack

by

 

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

Contact us