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]=volt_eq_creep(n,dom,kind)
function [x,solnum]=volt_eq_creep(n,dom,kind)
% Solves Volterra integral equations
% y(x)=int_a^x K(x,t)y(t)dt+f(x), x in [a,b]
% use: [x,solnum]=volt_eq_creep(32,0:0.05:0.2,2);
%

x=[];solnum=[]; k1 = 1000;k2 = 2000;eta = 10;d=length(dom);
for cont=1:d-1
    Vspec=volt(n,[dom(cont),dom(cont+1)],kind);x=[x;xv];
    F=x2t(myfun(xv),kind);
    for subcont=1:cont-1
        Fspec=fred(n,[dom(cont),dom(cont+1)],[dom(subcont),dom(subcont+1)],kind);
        F=F+Fspec*sol(:,subcont);
    end
    sol(:,cont)=(eye(n)-Vspec)\F;solnumnew=t2x(sol(:,cont),kind);
    solnum=[solnum;solnumnew];
end
myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Vspec,Vphys]=volt(n,domv,kind)
% Volterra operator with kernel K(x,t)
xv=pd(n,domv,kind);T=cpv(n,xv,domv);[~,J0]=prim(n,domv);
V=T*J0/T;
[YY,XX]=meshgrid(xv,xv);KW=mykernel(XX,YY);
Vphys=V.*KW;Vspec=T\Vphys*T;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Fspec,Fphys]=fred(n,domx,domt,kind)
% Fredholm operator with kernel K(x,t), x in domx, t in domt
xf=pd(n,domx,kind);[xt,w]=pd(n,domt,kind);T=cpv(n,xt,domt);
[TT,XX]=meshgrid(xt,xf);KW=mykernel(XX,TT);
Fphys=KW*diag(w);Fspec=T\Fphys*T;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res=mykernel(xk,tk)
    % x,t matrices
    % describes K(x,t), must be written by the user
    tauR = eta./k2;
    res = k2./tauR./(k1+k2).*exp(-(xk-tk)./tauR);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=myfun(xfun)
    % describes f(x) 
    % must be written by the user
    k1 = 1000;
    k2 = 2000;    
    f = ones(size(xfun))./(k1+k2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    % describes the output of the code
    % must be written by the user
    k1 = 1000;
    k2 = 2000;
    eta = 10;
    tauR = eta./k2;
    J = (1-k2/(k1+k2).*exp(-x./(1+k2/k1)./tauR))./k1;
    figure(1)
    plot(x,J,'-',x,solnum,'or');grid;xlabel('x');ylabel('y');
    legend('Analytical solution','Numerical solution');
    figure(2);semilogy(x,abs(J-solnum));grid;xlabel('x');ylabel('AbsErr');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us