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_newton_ex2(n,dom,kind,M,miu)
function [xg,solnumg]=ibvp_sys_newton_ex2(n,dom,kind,M,miu)
% Example:  y1'-y2=0, y2'+y1-miu y2=-miu y1^2 y2; x in [0,10]
%           y1(0)=2, y2(0)=0, miu=1 or miu=1000, (van der Pol equation)
% From: Matlab, vanderpoldemo.m
% Remark: uses Newton method with adaptive step
% call: [x,solnum]=ibvp_sys_newton_ex2(64,[0 10],2,16,1);
% call: [x,solnum]=ibvp_sys_newton_ex2(64,[0 3000],2,16,1000);
tic;xg=[];solnumg=[];solg=[];m=length(dom)-1;dyg=[];contg=0;step=1;
domloc=[dom(1),min(dom(1)+step,dom(2))];
myINIT;y0=yv(1,:);y=x2t(yv,kind);% starting approximation in physical space
while domloc(2)<dom(2)
    x=pd(n,domloc,kind);X=mult(n,domloc);J=prim(n,domloc);
    yv=[ones(n,1)*y0(1),ones(n,1)*y0(2)];y=x2t(yv,kind);
    dy=1;cont=1;
    while dy(cont) > 1.e-12,display(domloc(1));%display independent variable x
        if cont>30, break, end; 
        cont=cont+1;myDElin;myDEnonlin;myIC;cor=A\b;% this is the correction
        dy(cont)=norm(cor);if dy(cont)>100, break,end,
        y=reshape(cor,n,2)+y;yv=t2x(y,kind);
    end
    if cont>20 || norm(y(end-2:end,:))>1.e-12 || dy(cont)>100, step=step/2; 
        domloc=[domloc(1),domloc(1)+step]; 
    else  step=step*11/10;
    solnum = yv;
    solnumg=[solnumg;solnum];solg=[solg;y];
    xg=[xg;x];dyg=[dyg,dy];contg=contg+cont;
    y0=solnum(end,:);domloc=[domloc(2),min(domloc(2)+step,dom(2))];
    end;
end;toc;myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myINIT
    yv=[ones(n,1)*2,ones(n,1)*0];% initial approximation
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
    % must be written by the user
    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]*...
        ([zeros(n) -speye(n);speye(n),-miu*speye(n)]-DF);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
    % must be written by the user
    r=f(x,yv);F=x2t(r,kind);
    b=-(reshape(y,2*n,1)+[J,zeros(n);zeros(n),J]*...
        ([zeros(n) -speye(n);speye(n),-miu*speye(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)
    % must be written by the user
    rez=[zeros(n,1),-miu*yv(:,1).^2.*yv(:,2)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drez=df(x,yv)
    drez=[zeros(n,1),zeros(n,1);-2*miu*yv(:,1).*yv(:,2),-miu*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);semilogy(1:contg,dyg,'.');grid;
    title('History of iterations - corrections');
    figure(2);
    plot(xg,solnumg(:,1));grid;xlabel('x');ylabel('y_1(x)');
    title('Numerical solution');
    figure(3);
    plot(solnumg(:,1),solnumg(:,2));grid;xlabel('y_1(x)');ylabel('y_2(x)');
    title('Phase portrait of the solution');hold off
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us