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_ode_newton_test(n,dom,kind,M)
```function [x,solnum]=ibvp_ode_newton_test(n,dom,kind,M)
% Example: xy''+ 2y' = -xy^5, y'(0)=0, y(1)=sqrt(3)/2 (Emden's equation)
% Exact solution: solex = (1+x^2/3)^(-1/2);
% call: [x,solnum]=ibvp_ode_newton_test(32,[0 1],2,4);
%
tic;solnum=zeros(n,1);x=pd(n,dom,kind);X=mult(n,dom);D=deriv(n,dom);
yv=zeros(n,1);myINIT;y=x2t(yv,kind);% starting approximation in spectral space
dy=1;cont=1;
while dy(cont) > 1.e-14
if cont>100, break, end;
cont=cont+1;myDElin;myDEnonlin;cor=A\b;% this is the correction
dy(cont)=norm(cor);y=cor+y;yv=t2x(y,kind);
end;toc;solnum = yv;myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myINIT
% starting approximation in the physical space
yv=ones(n,1)*sqrt(3)/2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
% describes the linear part
% implements the initial/boundary conditions in the matrix A
dr=df(x,yv);DF=fact(x2t(dr,kind),M);A=X*D^2+2*D-DF;
T=cpv(n,[0,1],dom);A(n-1,:)=T(1,:)*D;A(n,:)=T(2,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
% describes the nonlinear part
% implements boundary conditions in the vector b
r=f(x,yv);b=-(X*D^2+2*D)*y+x2t(r,kind);b(n-1)=0;b(n)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rez=f(x,yv)
% describes the function f
rez=-x.*(yv.^5);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drez=df(x,yv)
% describes the Jacobian of the function f
drez=-5*x.*(yv.^4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
% describes the output of the code
figure(1);solex=1./sqrt(1+(x.^2)/3);
subplot(2,1,1);semilogy(abs(y),'.');grid;
title('Absolute value of the coefficients of the solution');
subplot(2,1,2);plot(x,solex-solnum);grid;
title('The error');xlabel('x');ylabel('err');
figure(2);
semilogy(1:cont,dy,'.');grid;
title('History of iterations - corrections');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end```