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:
 rfadd Addition of two rational functions error Display message and abort function. roots Find polynomial roots. unique Set unique.
This function is called by:
 mindist Finds minimal distance to the origin on the surface b'*x+x'*diag(g)*x=u

## 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