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.

Ex3( n)
function [ x,sol,dsol ] = Ex3( n)
% Example 3 
%
% u''=5./(1+exp(-8*(uu'-1))), 
% u(0)=int(u(s)ds,s=0..1)/2,
% u'(0)=int(u'(s)ds,s=0..1)/2=(u(1)-u(0))/2
%
% call: [ x,sol,dsol ] = Ex3( 64);
%

dom = [0,1];kind = 2;
[x,w] = pd(n,dom,kind);T = cpv(n,dom,dom);
D = deriv(n,dom);J=prim(n,dom);
options=optimset('Display','iter','TolFun',1.e-12,'MaxFunEvals',20000,...
    'MaxIter',1000,'Algorithm','Levenberg-Marquardt');
Uinit=x2t(1+x.^2,kind);

fact=5;U=fsolve(@sant,fact*Uinit,options);
sol(:,1)=t2x(U,kind);dsol(:,1)=t2x(D*U,kind);
subplot(2,2,1);plot(x,real(sol(:,1)),'-k',x,real(dsol(:,1)),'--k');grid;title('a');xlabel('t');ylabel('x,x''');
legend('x','x''');
T(1,:)*U-(T(2,:)-T(1,:))*J/2*U,T(1,:)*(D*U)-(T(2,:)-T(1,:))/2*U

fact=1;U=fsolve(@sant,fact*Uinit,options);
sol(:,1)=t2x(U,kind);dsol(:,1)=t2x(D*U,kind);
subplot(2,2,2);plot(x,real(sol(:,1)),'-k',x,real(dsol(:,1)),'--k');grid;title('b');xlabel('t');ylabel('y,y''');
legend('y','y''');
T(1,:)*U-(T(2,:)-T(1,:))*J/2*U,T(1,:)*(D*U)-(T(2,:)-T(1,:))/2*U

fact=0.1;U=fsolve(@sant,fact*Uinit,options);
sol(:,1)=t2x(U,kind);dsol(:,1)=t2x(D*U,kind);
subplot(2,2,3);plot(x,real(sol(:,1)),'-k',x,real(dsol(:,1)),'--k');grid;title('c');xlabel('t');ylabel('z,z''');
legend('z','z''');
T(1,:)*U-(T(2,:)-T(1,:))*J/2*U,T(1,:)*(D*U)-(T(2,:)-T(1,:))/2*U

u=linspace(0,12.5,1000);
F=5./u./(1+exp(-8*(u.^2/2-1)))-2;
G=5./u./(1+exp(-8*(3*u.^2/2-1)))-2/3;
H=5./u./(1+exp(-8*(16*u-1)))-2/3;
subplot(2,2,4);plot(u,F,'-k',u,G,'--k',u,H,':k');axis([0 12.5 -2.2 5]);grid;
legend('F','G','H');xlabel('x');ylabel('F,G,H');title('d');

    function F=sant(U)
        % U coefficients
        u=t2x(U,kind);du=t2x(D*U,kind);
        f=x2t(5./(1+exp(-8*(u.*du-1))),kind);
        A=D^2;
        A(n-1,:)=T(1,:)-(T(2,:)-T(1,:))*J/2;f(n-1)=0;
        A(n,:)=T(1,:)*D-(T(2,:)-T(1,:))/2;f(n)=0;
        F=A*U-f;
    end
end

Contact us