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.

[x,solnum]=ibvp_int_succapprox_ex1(n,dom,kind)
function [x,solnum]=ibvp_int_succapprox_ex1(n,dom,kind)
% Example: y''(x)+xy'(x)-xy(x)=-x-2x^2+int_0^1{6xy^2(t)}dt+
%                +int_0^x{2y(t)}dt, x in [0,1]
%          y(0)=0, y'(0)=1
% From: M. Dadkhah, M. T. Kajani and S. Mahdavi, Proceedings of the 6th IMT-GT 
%       Conference on Mathematics, Statistics and its Applications (ICMSA2010)
%       Universiti Tunku Abdul Rahman, Kuala Lumpur,Malaysia, pp. 738-744
% Exact solution: solex=x
% call: [x,solnum]=ibvp_int_succapprox_ex1(32,[0 1],2); 
%
tic;
[x,w]=pd(n,dom,kind);D=deriv(n,dom);X=mult(n,dom);T=cpv(n,x,dom);
[YY,XX]=meshgrid(x,x);myINIT;dy=1;cont=1;KF=zeros(n,1);KV=KF;
[~,Fphys]=fred(n,dom,kind);[~,Vphys]=volt(n,dom,kind);
while dy(cont) > 1.e-12
    if cont>100, break, end; 
    cont=cont+1; yoval=T*yo;
    KV=Vphys*yoval;% values of the Volterra integral
    KF=Fphys*(yoval.^2);% values of the Fredholm integral   
    myDE;myBC;yn=A\b;dy(cont)=norm(yn-yo);yo=yn;
end;
toc;sol = yn;solnum=t2x(sol,kind);myOUT;
function [Fspec,Fphys]=fred(~,~,~)
% Fredholm operator with kernel K(x,t)
KWF=myfredkernel(XX,YY);
Fphys=KWF*diag(w);Fspec=T\Fphys*T;
end
function [Vspec,Vphys]=volt(n,dom,~)
% Volterra operator with kernel K(x,t)
[~,J0]=prim(n,dom);
V=T*J0/T;
KWV=myvoltkernel(XX,YY);
Vphys=V.*KWV;Vspec=T\Vphys*T;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res=myfredkernel(x,~)
    % x,t matrices
    res=6*x;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res=myvoltkernel(~,~)
    % x,t matrices
    res=2*ones(n);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
    A=D^2+X*D-X;f=-x-2*x.^2+KF+KV;b=x2t(f,kind);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
    Tbc=cpv(n,dom(1),dom);A(n-1,:)=Tbc;b(n-1)=0;A(n,:)=Tbc*D;b(n)=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myINIT
    % describes the initial approximation in spectral form
    yo=x2t(sin(x),kind);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    solex=x;
    figure(1);
    subplot(2,1,1);semilogy(abs(sol),'.');grid;
    title('Absolute value of the coefficients of the solution');
    subplot(2,1,2);plot(x,solnum);grid;xlabel('x');ylabel('y(x)');
    title('The numerical solution');
    figure(2);
    semilogy(1:cont,dy,'.');grid;title('History of iterations - corrections');
    figure(3);plot(x,abs(solnum-solex));grid;xlabel('x');ylabel('err'); 
    title('The error');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


Contact us