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.

[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