Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

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```