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.

[time,x,y]=ex6(k,v0x,v0y)
function [time,x,y]=ex6(k,v0x,v0y)
% Shampine's Example ch2ex3.m from 
% http://www.radford.edu/~thompson/sodeswm/mfiles/mfiles.html or
% Solving ODEs with Matlab, Cambridge University Press, 2003
% L.F. Shampine, I. Gladwell, S. Thompson, pp. 98-102, Example 2.3.5
% use [time,x,y]=ex6(0.32,0,0);
%
fill([0 0 1],[0 1 0],[0.8 0.8 0.8]);
axis([-0.1 1.1 0 4]);
hold on
plot(0,2,'ro');
n=64;dom=[0 1];kind=2;g=9.81;tol=1.e-13;
b=zeros(n,2);b(1,2)=-g;
x0=0;xp0=v0x;y0=2;yp0=v0y;%viteza initiala
time=[];x=[];y=[];
%tic;
for j=1:50
    D=deriv(n,dom);T=cpv(n,dom(1),dom);t=pd(n,dom,kind);
    A=D^2;A(n-1,:)=T;A(n,:)=T*D;b(n-1:n,:)=[x0,y0;xp0,yp0];
    sol=A\b;solnum=t2x(sol,kind);dsol=D*sol;dsolnum=t2x(dsol,kind);
    J=find(solnum(:,1)<=1);
% find the touch time tk
    c=sol(:,1)+sol(:,2)-[2;zeros(n-1,1)];
    [z,fz,m]=chebroots(c,dom,dom,tol);
    z=real(z);z=z(abs(fz)<tol);I=length(z);
    if I==0
        plot(solnum(J,1),solnum(J,2),'r');hold on;%pause(0.2);
        x0=solnum(end,1);xp0=dsolnum(end,1);
        y0=solnum(end,2);yp0=dsolnum(end,2);
    else tk=z(1);
        Tr=cpv(n,tk,dom);
        x0=Tr*sol(:,1);xp0=-k*Tr*(D*sol(:,2));
        y0=Tr*sol(:,2);yp0=k*Tr*(D*sol(:,1));
        J1=find(t<tk);J2=intersect(J,J1);
        if x0<=1,
            plot([solnum(J2,1);x0],[solnum(J2,2);y0],'r',x0,y0,'og');
            hold on;%pause(0.2);
        else 
            plot(solnum(J2,1),solnum(J2,2),'r');hold on;%pause(0.2);
        end
        if x0>=1, break;end
        if ~isempty(time)&&(tk-time(end))/time(end)<0.001, 
            display('Cluster, end bouncing');break;
        end
        time=[time;tk];x=[x;x0];y=[y;y0];dom=[tk,tk+1];
    end
end
%toc;
grid;xlabel('x');ylabel('y');hold off;
display([time,x,y]);

Contact us