Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.fr>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Common tangent
Date: Wed, 28 Nov 2007 17:07:45 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 128
Message-ID: <fik791$1f9$1@fred.mathworks.com>
References: <fihu0m$quf$1@fred.mathworks.com> <fiivds$lqs$1@fred.mathworks.com> <fij6n6$9h7$1@fred.mathworks.com> <fijhn9$bfi$1@fred.mathworks.com> <fik01e$csg$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.fr>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1196269665 1513 172.30.248.37 (28 Nov 2007 17:07:45 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 28 Nov 2007 17:07:45 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:439830


A little program for fun. Zeros solver is dirty written, one
can do much better.

Bruno

% Input function
x=[-5:0.5:5];
f=sin(x) + 0.01*x.^3;
g=sin(x/2+1) - 2;

%
% Call function
%
commontangent(x, f, g);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function commontangent(x, f, g)
%
% Solve and graphic output for common tangent
%

% Spline
F=spline(x,f);
G=spline(x,g);
% Derivative
DF=splinederive(F);
DG=splinederive(G);
% Constant term
AF=coef0(F,DF);
AG=coef0(G,DG);

% fine grid
xx=linspace(min(x),max(x),100);

af=ppval(AF,xx); bf=ppval(DF,xx);
ag=ppval(AG,xx); bg=ppval(DG,xx);

cf=[af./sqrt(1+bf.^2); bf; xx];
cg=[ag./sqrt(1+bg.^2); bg; xx];

%
% Solve for common tangent; dirty way
%
c=poly2poly(cf, cg);

%
% Graphic output
%
figure(1);
clf;

subplot(2,1,1);
hf=plot(x,f,'b+',...
        xx,ppval(F,xx),'b',...
        x,g,'g+',...
        xx,ppval(G,xx),'g');
hold on;

for n=1:size(c,2)
    cn=c(:,n);
    b=cn(2);
    a=cn(1)*sqrt(1+b^2);
    x2=cn([3 6]); 
    ht=plot(x2, polyval([b a], x2), 'ko-');
end
h=[hf(1) hf(3) ht(1)];
legend(h, 'f', 'g', 'common tangents');

subplot(2,1,2);
h=plot(cf(1,:),cf(2,:),'b',...
       cg(1,:),cg(2,:),'g',...
       c(1,:), c(2,:), 'k.');
legend(h, 'f-tangent', 'g-tangent', 'common tangents'); 

end

% derivative of a spline
function d=splinederive(pp)
d=pp;
d.order=pp.order-1;
order=d.order:-1:1;
d.coefs=pp.coefs(:,1:d.order).*repmat(order,pp.pieces,1);
end

% constant coefficient of the tangent
function A=coef0(F,DF)
A=F;
z=zeros(DF.pieces,1);
P=[DF.coefs z] + ...
  [z repmat(DF.breaks(1:end-1)',1,DF.order).*DF.coefs];
A.coefs=F.coefs-P;
end

% intersection of two line segments
function c=seg2seg(s1, s2)

M=[s1(1:2,2)-s1(1:2,1) -s2(1:2,2)+s2(1:2,1)];
rhs=s2(1:2,1)-s1(1:2,1);
sol=M\rhs;
if all(sol>=0 & sol<=1)
    c=[s1*[1-sol(1) sol(1)]'; ...
       s2*[1-sol(2) sol(2)]']; 
else
    c=[];
end

end

% intersection of a q segment and a polygone
function c=seg2poly(s1, P)
c=[];
for i=2:size(P,2)
    c=[c seg2seg(s1, P(:,i-1:i))]; %#ok
end
end

% intersection of two polygones
function c=poly2poly(P1, P2)
c=[];
for i=2:size(P1,2)
    c=[c seg2poly(P1(:,i-1:i), P2)];%#ok
end
end