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]=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

```