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_test(n,dom,kind,y0)
function [xg,solnumg]=ibvp_sys_succapprox_test(n,dom,kind,y0)
% Example: y1'=Ay1-B y1 y2;y2'= -Cy2+B y1 y2; x in [0,100]
%          y1(0)=0.5, y2(0)=0.5, A=0.1, B=0.2, C=0.05
%          (Lotka-Volterra system)
% test: invariant=By1+By2-Clog y1-Alog y2
% call: [x,solnum]=ibvp_sys_succapprox_test(32,[0:10:100],2,[0.5,0.5]);
%
tic;xg=[];solnumg=[];solg=[];p=length(dom)-1;dyg=[];contg=0;m=length(y0);
A=spalloc(p*m,p*m,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)-0.1*J,zeros(n);zeros(n),speye(n)+0.05*J];
    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).*yv(:,2);
    b=[J,zeros(n);zeros(n),J]*reshape(x2t([-0.2*work,0.2*work],kind),2*n,1);
    b(1)=y0(1);b(n+1)=y0(2); 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    % describes the output of the code
    invar=0.2*solnumg(:,1)+0.2*solnumg(:,2)-0.05*log(solnumg(:,1))-...
        0.1*log(solnumg(:,2));
    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(xg,invar);xlabel('xg');ylabel('The invariant');
    title('The invariant of the solution');
    figure(4);
    xx=linspace(dom(1),dom(end),100000);
    f1 = barycheb(xx,solnumg(:,1),xg,kind);f2 = barycheb(xx,solnumg(:,2),xg,kind);
    plot(f1,f2);xlabel('y1(x)');ylabel('y2(x)');grid;
    title('The phase portrait of the solution');
    figure(5);spy(A);title('spy(A)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us