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.

[xx,solnum]=ibvp_ode_split_ex3(n,dom,kind)
function [xx,solnum]=ibvp_ode_split_ex3(n,dom,kind)
% Example: -ep y'''+4 y'- y = 0.7-(x>0.5)*1.3, x in [0,1]
%          y(0)=1, y'(0)=0, y'(1)=0, ep=1/2^8
% From: T. Valanarasu, N. Ramanujam, 
%       Novi Sad J. Math., Vol. 37, No. 2, 2007, 41-57) 
% call: [x,solnum]=ibvp_ode_split_ex3(32,[0 0.5 1],2);
% or    [x,solnum]=ibvp_ode_split_ex3(32,[0 0.5 1],1);
%           where we have no gridpoints at 0.5
%
tic;ep=1/2^8;
m=length(dom);
A=spalloc(n*(m-1),n*(m-1),100*n*(m-1));b=zeros(n*(m-1),1);xx=[];
xs=pd(n,[-1 1],kind);Js=prim(n,[-1,1]);Ds=deriv(n,[-1,1]);
for k=1:m-1
    pdom=[dom(k),dom(k+1)];l=(pdom(2)-pdom(1))/2;med=(pdom(2)+pdom(1))/2;
    x=l*xs+med;J=l*Js;xx=[xx,x];
    myDE;cff(:,k)=cf;
    E=zeros(m-1);E(k,k)=1;AA(1:2,:)=0;A=A+kron(E,AA);
    EE=zeros(m-1,1);EE(k)=1;bb(1:2)=0;b=b+kron(EE,bb);
end;
myBC; sol=A\b;% the solution in spectral form
toc;solnum=[];dsolnum=[];
for i=1:m-1
    solnum = [solnum,t2x(sol((i-1)*n+1:i*n),kind)];
    dsolnum = [dsolnum,t2x(2*Ds/(dom(i+1)-dom(i))*sol((i-1)*n+1:i*n),kind)];
end
% the solution in physical form
x=reshape(xx,(m-1)*n,1);
solnumv=reshape(solnum,(m-1)*n,1);dsolnumv=reshape(dsolnum,(m-1)*n,1);
myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
    if x(n-1)<0.5, x(n)=x(n)-eps; else x(1)=x(1)+eps;end
    f=0.7-(x>0.5)*1.3;cf=x2t(f,kind);
    AA=-ep*speye(n)+4*J^2-J^3;bb=J^3*cf;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
    T=cpv(n,[dom(1),dom(m)],[dom(1),dom(m)]);% T(1,:) for xc=0, T(2,:) for xc=1
    Tl=T(1,:);Tr=T(2,:);
    A(1,1:n)=Tl;A(2,1:n)=Tl*2*Ds/(dom(2)-dom(1));
    A(3,(m-2)*n+1:(m-1)*n)=Tr*2*Ds/(dom(m)-dom(m-1));
    b(1)=1;b(2)=0;b(3)=0;
    for j=2:m-1 
        fl=(dom(j)-dom(j-1))/2;fr=(dom(j+1)-dom(j))/2;
        A((j-1)*n+1,(j-2)*n+1:(j-1)*n)=Tr;
        A((j-1)*n+1,(j-1)*n+1:j*n)=-Tl;b((j-1)*n+1)=0;
        A((j-1)*n+2,(j-2)*n+1:(j-1)*n)=Tr*Ds/fl;
        A((j-1)*n+2,(j-1)*n+1:j*n)=-Tl*Ds/fr;b((j-1)*n+2)=0;
        A((j-1)*n+3,(j-2)*n+1:(j-1)*n)=Tr*(Ds/fl)^2;
        A((j-1)*n+3,(j-1)*n+1:j*n)=-Tl*(Ds/fr)^2;b((j-1)*n+3)=0;
    end    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    figure(1);
    subplot(2,1,1);semilogy(abs(sol),'.');grid;
    title('The absolute value of the coefficients of the solution');
    subplot(2,1,2);plot(reshape(cff,(m-1)*n,1),'.');grid;
    title('The coefficients of the rhs');
    figure(2);
    subplot(2,1,1);plot(x,solnumv,'b',dom,zeros(length(dom),1),'.r');
    grid;xlabel('x');ylabel('y(x)');title('The numerical solution');
    subplot(2,1,2);plot(x,dsolnumv,'b',dom,zeros(length(dom),1),'.r');
    grid;xlabel('x');ylabel('y''(x)');
    title('The derivative of the numerical solution');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us