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.

ibvp_ode_newton_ex3.m
function [x,solnum,A]=ibvp_ode_newton_ex3(n,dom,kind,M)
% Example: y''' = - (y+y^2)^(1/3), x in [0,1]
%          y(0)-y'(0)=0,y'(3/10)=0, y(1)=0
% From: A. P. PALAMIDES, N. M. STAVRAKAKIS,
% Electronic Journal of Differential Equations, 2010(155), pp. 112.
% call: [x,solnum]=ibvp_ode_newton_ex3(32,[0 1],1,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-13
    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
    yv=(-x.^3-23/10*x.^2+33/20*x+33/20);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
    dr=df(x,yv);DF=fact(x2t(dr,kind),M);A=D^3-DF;
    T=cpv(n,[0,3/10,1],dom);
    A(n-2,:)=T(1,:)-T(1,:)*D;A(n-1,:)=T(2,:)*D;A(n,:)=T(3,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
    r=f(x,yv);b=-D^3*y+x2t(r,kind);b(n-2)=0;b(n-1)=0;b(n)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rez=f(x,yv)
    rez=-nthroot((yv+yv.^2),3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drez=df(x,yv)
    drez=-1/3*(1+2*yv)./nthroot((yv+yv.^2).^2,3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    figure(1);
    subplot(2,1,1);semilogy(abs(y),'.');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');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


Contact us