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.

[xg,solnumg]=ibvp_sys_succapprox_ex1(n,dom,kind,y0)
function [xg,solnumg]=ibvp_sys_succapprox_ex1(n,dom,kind,y0)
% 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
%          (the Brusselator)
% test: fixed point at (B,A/B)
% call: [x,solnum]=ibvp_sys_succapprox_ex1(64,[0,1,2:2:200],2,[1.5,3]);
%
tic;xg=[];solnumg=[];solg=[];p=length(dom)-1;dyg=[];contg=0;m=length(y0);
A=spalloc(m*n,m*n,10*p*m);
for i=1:p
    domloc=[dom(i),dom(i+1)];
    x=pd(n,domloc,kind);X=mult(n,domloc);J=prim(n,domloc);xg=[xg;x];
    myDElin;y=zeros(m*n,1);cont=1;dy=1;
    for k=1:m,y((k-1)*n+1)=2*y0(k);end
    while dy(cont) > 1.e-12
        if cont>100, break, end; 
        cont=cont+1;myDEnonlin;sol=A\b;dy(cont)=norm(sol-y);y=sol;
    end
    solnum = t2x(reshape(sol,n,m),kind);
    solnumg=[solnumg;solnum];solg=[solg;sol];dyg=[dyg,dy];contg=contg+cont;
    y0=solnum(end,:);
end;toc;myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
    % describes the linear part of the differential system
    % implements the initial conditions in the matrix A
    A=[speye(n)+2.9*J,zeros(n);-1.9*J,speye(n)];
    T=cpv(n,domloc(1),domloc);A(1,:)=[T,zeros(1,n)];A(n+1,:)=[zeros(1,n),T];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
    % describes the nonlinear part of the differential system
    % implements the initial conditions in the vector b
    yv=t2x(reshape(y,n,2),kind);work=(yv(:,1).^2).*yv(:,2);
    b=[J,zeros(n);zeros(n),J]*reshape(x2t([1+work,-work],kind),2*n,1);
    b(1)=y0(1);b(n+1)=y0(2); 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    % describes the output of the code
    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('Phase portrait');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us