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_int_ex4.m
function [x,solnum]=ibvp_ode_int_ex4(n,dom,kind)
% Example:  x^4y^(4)- 4x^3y'''+x^2(2-x^2)y''+2x(x^2-12)y'+2(12-x^2)y=2x^5, x in [1,11]                                                                                       
%           y(1)= 1+e+1/e, y'(1)=2e, y(11)=-1199+11e^11+11e^(-11),
%           y'(11)=-340+12e^11-10e^(-11)
% From: N. Mai-Duy, International Journal for Numerical Methods 
%       in Engineering, 62, 6, 2005, pp. 824852
% Exact solution: solex=x+x^2-x^3+xexp(x)+xexp(-x)
% call: [x,solnum]=ibvp_ode_int_ex4(64,[1 11],2);
% 
tic;
A=zeros(n);b=zeros(n,1);
x=pd(n,dom,kind);X=mult(n,dom);J=prim(n,dom);D=deriv(n,dom);
myDE;myBC;sol=A\b;% the solution in spectral form
toc;
solnum = t2x(sol,kind);% the solution in physical form
myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
    f=2*x.^5;
    A4=X^4-16*J*X^3+72*J^2*X^2-96*J^3*X+24*J^4;
    A3=J*(-4*X^3+36*J*X^2-72*J^2*X+24*J^3);
    A2=J^2*(12*X^2-X^4-2*J*(24*X-4*X^3)+J^2*(24*speye(n)-12*X^2));
    A1=J^3*(2*X^3-24*X-J*(6*X^2-24*speye(n)));A0=J^4*(24*speye(n)-2*X^2);
    A=A4+A3+A2+A1+A0;b=J^4*x2t(f,kind);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
    T=cpv(n,dom,dom);% T(1,:) for xc=0, T(2,:) for xc=1
    A(1,:)=T(1,:);b(1)=1+exp(1)+exp(-1);
    A(2,:)=T(1,:)*D;b(2)=2*exp(1);
    A(3,:)=T(2,:);b(3)=-1199+11*exp(11)+11*exp(-11);
    A(4,:)=T(2,:)*D;b(4)=-340+12*exp(11)-10*exp(-11);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
    solex=x+x.^2-x.^3+x.*exp(x)+x.*exp(-x);
    figure(1)
    subplot(2,1,1);semilogy(abs(sol),'.');grid;
    title('The absolute value of the coefficients of the solution');
    subplot(2,1,2);plot(x,solnum);grid;xlabel('x');ylabel('y(x)');
    title(['The solution (error = ',num2str(norm(solnum-solex)),')']);
    figure(2);semilogy(x,abs(solnum-solex));title('The errors');grid;
    xlabel('x');ylabel('err');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


Contact us