Home > wafo > trgauss > private > mindist3.m

mindist3

PURPOSE ^

Finds point of minimal distance to the origin on the surface b'*x+g'*x.^2=u.

SYNOPSIS ^

[x0,pmx]=mindist3(g,b,u)

DESCRIPTION ^

  MINDIST3 Finds point of minimal distance to the origin on the surface b'*x+g'*x.^2=u. 
  
    CALL:   [x0,pmx] = mindist3(g,b,u); 
  
  Returns a point of minimal distance to the origin on the surface 
  b'*x+g'*x.^2=u

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [x0,pmx]=mindist3(g,b,u) 
002 % MINDIST3 Finds point of minimal distance to the origin on the surface b'*x+g'*x.^2=u. 
003 % 
004 %   CALL:   [x0,pmx] = mindist3(g,b,u); 
005 % 
006 % Returns a point of minimal distance to the origin on the surface 
007 % b'*x+g'*x.^2=u 
008  
009 g=g(:); 
010 b=b(:); 
011 d0=length(g); 
012 if length(b)~=d0 
013    error('g and b should have the same length') 
014 end 
015 % If u==0, x0=0 is the unique solution 
016 if u==0 
017    x0=zeros(d0,1); 
018    pmx=zeros(d0,1); 
019    return 
020 end 
021  
022 nulldim=(g==0)&(b==0);% If b(i)=g(i)=0, x0(i)=0 
023  
024 if sum(nulldim)==d0% If all b's and g's are zeros and u=0, only zeros are returned 
025    if u==0 
026       x0=zeros(d0,1); 
027       pmx=zeros(d0,1); 
028       return 
029    else 
030       error('u is out of range') 
031    end 
032 end 
033  
034 g(nulldim)=[]; 
035 b(nulldim)=[]; 
036 d=length(b); 
037 if (all(g>0)|all(g<0))&u==-.25*sum(b.^2./g) 
038    x00=-.5*b./g; 
039    pmx=zeros(d,1); 
040    x0=zeros(d0,1); 
041    x0(~nulldim)=x00; 
042 end 
043 if all(g>0)&u<-.25*sum(b.^2./g) 
044    error('u out of range') 
045 end 
046 if all(g<0)&u>-.25*sum(b.^2./g) 
047    error('u out of range') 
048 end 
049 pmx0=zeros(d,1); 
050 x00=[inf;zeros(d-1,1)]; 
051 I0=b==0; 
052 I1=b~=0; 
053 gI0=g(I0); 
054 gI0unique=unique(gI0); 
055 for k=1:length(gI0unique) 
056    if all(g(I1)-gI0unique(k)) 
057       xp=zeros(d,1); 
058       xp(I1)=.5*b(I1)./(gI0unique(k)-g(I1)); 
059       up=u-sum(b.*xp+g.*xp.^2,1); 
060       if sign(up)*sign(gI0unique(k))>=0 
061      xp(min(find(g==gI0unique(k))))=sqrt(up/gI0unique(k)); 
062      if sum(xp.^2)<sum(x00.^2) 
063         x00=xp; 
064         pmx0=zeros(d,1); 
065         pmx0(g==gI0unique(k))=1; 
066      end  
067       end 
068    end  
069 end 
070 % Construction of the rational function: 
071 bet=b(I1); 
072 gam=g(I1); 
073 if sum(I1)>0 
074    R=[bet(1)^2*[0 2 -gam(1)];1 -2*gam(1) gam(1)^2]; 
075    for k=2:sum(I1); 
076       R=rfadd(R,[bet(k)^2*[0 2 -gam(k)];1 -2*gam(k) gam(k)^2]); 
077    end 
078    lam=roots(R(1,:)-4*u*R(2,:)); 
079    lam=real(lam(abs(real(lam))>(1e6)*abs(imag(lam)))); 
080    x1=zeros(d,1); 
081    for k=1:length(lam) 
082       x1(I1)=.5*bet./(lam(k)-gam); 
083       if sum(x1.^2)<sum(x00.^2); 
084      x00=x1; 
085      pmx0=zeros(size(pmx0)); 
086       end 
087    end 
088 end 
089  
090  
091 x0=zeros(d0,1); 
092 x0(~nulldim)=x00; 
093 pmx=zeros(d0,1); 
094 pmx(~nulldim)=pmx0; 
095

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003