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.

[x,solnum]=ibvp_sys_newton_ex1(n,dom,kind,M)
function [x,solnum]=ibvp_sys_newton_ex1(n,dom,kind,M)
% Example: y1'=B-(A+1)y1+y1^2 y2;y2'=Ay1-y1^2 y2; x in [0,200]
%          y1(0)=1.5, y2(0)=3, A=1.9, B=1 (Brusselator system) 
% Test: fixed point at (B,A/B)
% call: [x,solnum]=ibvp_sys_newton_ex1(64,[0 1 2:2:200],2,16);
% 
tic;coefA=1.9;coefB=1;
xg=[];solnumg=[];solg=[];m=length(dom)-1;dyg=[];contg=0;
myINIT;y0=yv(1,:);% starting approximation in physical space
h = waitbar(0,'Please wait...');
for i=1:m
    waitbar(i/m,h);
    domloc=[dom(i),dom(i+1)];
    x=pd(n,domloc,kind);X=mult(n,domloc);J=prim(n,domloc);xg=[xg;x];
    yv=[ones(n,1)*y0(1),ones(n,1)*y0(2)];y=x2t(yv,kind);% starting approximation
    dy=1;cont=1;
    while dy(cont) > 1.e-12
        if cont>100, break, end; 
        cont=cont+1;myDElin;myDEnonlin;myIC;cor=A\b;% this is the correction
        dy(cont)=norm(cor);y=reshape(cor,n,2)+y;yv=t2x(y,kind);
    end
solnum = yv;solnumg=[solnumg;solnum];solg=[solg;y];dyg=[dyg,dy];
contg=contg+cont;y0=solnum(end,:);
end;toc;
close(h);myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myINIT
    yv=[ones(n,1)*1.5,ones(n,1)*3];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
    dr=df(x,yv);
      DF=[fact(x2t(dr(1:n,1),kind),M),fact(x2t(dr(1:n,2),kind),M);...
          fact(x2t(dr(n+1:2*n,1),kind),M),fact(x2t(dr(n+1:2*n,2),kind),M)];
    A=[speye(n) zeros(n);zeros(n) speye(n)]+[J,zeros(n);zeros(n),J]*...
        ([(coefA+1)*speye(n) zeros(n);-coefA*speye(n),zeros(n)]-DF);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
    r=f(x,yv);F=x2t(r,kind);
    b=-(reshape(y,2*n,1)+[J,zeros(n);zeros(n),J]*...
        ([(coefA+1)*speye(n) zeros(n);-coefA*speye(n),zeros(n)]*reshape(y,2*n,1)-...
        reshape(F,2*n,1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myIC
    T=cpv(n,dom(1),dom);A(1,:)=0;A(1,1:n)=T;A(n+1,:)=0;A(n+1,n+1:2*n)=T;
    b(1)=0;b(n+1)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rez=f(x,yv)
    rez=[coefB+yv(:,1).^2.*yv(:,2),-yv(:,1).^2.*yv(:,2)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drez=df(x,yv)
    drez=[2*yv(:,1).*yv(:,2),yv(:,1).^2;-2*yv(:,1).*yv(:,2),-yv(:,1).^2];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    figure(1);
    subplot(2,1,1);semilogy(abs(solg),'.');
    title('Absolute value of the coefficients of y_1, y_2');grid;
    subplot(2,1,2);plot(xg,solnumg(:,1),'b',xg,solnumg(:,2),'r');
    title('Numerical solution');legend('y1(x)','y2(x)');
    xlabel('x');ylabel('y1,2(x)');grid;
    figure(2);
    semilogy(1:contg,dyg,'.');grid;title('History of iterations - corrections');
    figure(3);
    plot(solnumg(:,1),solnumg(:,2),'b',1,1.9,'.r');grid;
    xlabel('y1(x)');ylabel('y2(x)');title('The phase portrait');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


Contact us