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.

ibvp_sys_ex2x2.m
function [x,solnum]=ibvp_sys_ex2x2(n,dom,kind,y0)
% Example: y1''+0.3y1'+8y1-4y2=0;y2''+0.3y2'+4y2-4y1=0, x in [0,20]
%          y1(0)=-1, y2(0)=1, y1'(0)=0, y2'(0)=0
% From: T. A. Driscoll, F. Bornemann, L. N. Trefethen, 
%       BIT Numerical Mathematics, (2008) 48: 701723
% call: [x,solnum]=ibvp_sys_ex2x2(64,[0,20],2,[-1,1,0,0]);
%
tic;m=length(y0)/2;A=spalloc(m*n,m*n,5*m*n);b=zeros(m*n,1);
x=pd(n,dom,kind);D=deriv(n,dom);
myDE;myBC;sol=reshape(A\b,n,m);% the solution in spectral space
solnum = t2x(sol,kind);dsolnum=t2x(D*sol,kind);% the solution in physical space
toc;myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
    b=zeros(m*n,1);
    A=[D^2+0.3*D+8*speye(n),-4*speye(n);-4*speye(n),D^2+0.3*D+4*speye(n)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
    T=cpv(n,dom(1),dom);
    A(n-1,:)=0;A(n-1,1:n)=T;b(n-1)=y0(1);
    A(n,:)=0;A(n,1:n)=T*D;b(n)=y0(3);
    A(2*n-1,:)=0;A(2*n-1,n+1:2*n)=T;b(2*n-1)=y0(2);
    A(2*n,:)=0;A(2*n,n+1:2*n)=T*D;b(2*n)=y0(4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    xx=linspace(dom(1),dom(2),10000);
    fx1 = barycheb(xx,solnum(:,1),x,kind);
    fx2 = barycheb(xx,solnum(:,2),x,kind);
    subplot(2,1,1);semilogy(abs(sol),'.');
    title('Absolute value of the coefficients of y_1 and y_2');grid;
    subplot(2,1,2);plot(xx,fx1,'b',xx,fx2,'g');title('Numerical solution');
    xlabel('x');ylabel('y_{1,2}(x)');grid;legend('y_1(x)','y_2(x)');
    energy=(dsolnum(:,1).^2+dsolnum(:,2).^2)/2+...
        4*(solnum(:,1).^2+(solnum(:,2)-solnum(:,1)).^2)/2;
    display(energy(end));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us