{VERSION 3 0 "IBM INTEL NT" "3.0" }
{USTYLETAB {CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 
0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }
{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "2
D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 256 "
" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 
0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 
"" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 
0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 
0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 
{CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 
0 0 0 0 -1 0 }{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 
0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 257 44 "Appendix A
--A Step in the Dual Newton Method" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
{TEXT 258 26 "Aleksandar Donev, 12/17/00" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 127 "This Maple worksheet illustrates the basic step in the D
ual Newton Method for convex network optimization when the arcs have a
 " }{TEXT 267 33 "superconductor-like cost function" }{TEXT -1 3 ".  \+
" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "restart: " "6#%(restartG
" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "with(linalg):" "6#-%%with
G6#%'linalgG" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "with(plots):
" "6#-%%withG6#%&plotsG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }
{TEXT 256 52 "Part I: Making the network node-arc incidence matrix" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "The actual node-arc incidence mat
rix of this graph is composed of two parts:" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 48 "The basis part corresponding to a spanning tree:" }}}
{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "A[B] := matrix([[1, 0, -1, 0, \+
1, 0], [0, 0, 0, 1, -1, -1], [0, 0, 0, -1, 0, 0], [0, 0, 0, 0, 0, 1], \+
[0, 0, 1, 0, 0, 0], [-1, -1, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]]):" "6#>&
%\"AG6#%\"BG-%'matrixG6#7)7(\"\"\"\"\"!,$\"\"\"!\"\"F.\"\"\"F.7(F.F.F.
\"\"\",$\"\"\"F1,$\"\"\"F17(F.F.F.,$\"\"\"F1F.F.7(F.F.F.F.F.\"\"\"7(F.
F.\"\"\"F.F.F.7(,$\"\"\"F1,$\"\"\"F1F.F.F.F.7(F.\"\"\"F.F.F.F." }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 52 "A non-basis part corresponding to \+
the non-tree arcs:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "A[N] :=
 matrix([[0, 0, 0], [0, 0, 0], [0, 1, 0], [0, -1, 0], [1, 0, -1], [0, \+
0, 1], [-1, 0, 0]]):" "6#>&%\"AG6#%\"NG-%'matrixG6#7)7%\"\"!F-F-7%F-F-
F-7%F-\"\"\"F-7%F-,$\"\"\"!\"\"F-7%\"\"\"F-,$\"\"\"F47%F-F-\"\"\"7%,$
\"\"\"F4F-F-" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 52 "And the full matr
ix is the concatenation of the two:" }}}{EXCHG {PARA 0 "> " 0 "" 
{XPPEDIT 19 1 "A[G]:=concat(A[B],A[N]):" "6#>&%\"AG6#%\"GG-%'concatG6$
&F%6#%\"BG&F%6#%\"NG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 113 "This mat
rix has rank 6 even though it has 7 rows. So we should throw away one \+
row corresponding to the root node:" }}}{EXCHG {PARA 0 "> " 0 "" 
{XPPEDIT 19 1 "rank(A[G]);" "6#-%%rankG6#&%\"AG6#%\"GG" }}{PARA 11 "" 
1 "" {XPPMATH 20 "6#\"\"'" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "Take
 Root to be 1:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "Root := 1:
" "6#>%%RootG\"\"\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 50 "So the act
ual node-arc incidence matrix we use is:" }}}{EXCHG {PARA 0 "> " 0 "" 
{XPPEDIT 19 1 "A:=delrows(A[G],Root..Root);" "6#>%\"AG-%(delrowsG6$&F$
6#%\"GG;%%RootGF," }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"AG-%'matrixG6
#7(7+\"\"!F*F*\"\"\"!\"\"F,F*F*F*7+F*F*F*F,F*F*F*F+F*7+F*F*F*F*F*F+F*F
,F*7+F*F*F+F*F*F*F+F*F,7+F,F,F*F*F*F*F*F*F+7+F*F+F*F*F*F*F,F*F*" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "This is now indeed a full-rank mat
rix:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "rank(A);" "6#-%%rankG
6#%\"AG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"'" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 37 "The number of arcs and nodes is thus:" }}}{EXCHG 
{PARA 0 "> " 0 "" {XPPEDIT 19 1 "m:=9: n:=6:" "6#C$>%\"mG\"\"*>%\"nG\"
\"'" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 257 54 "Part II) C
ost and conjugate functions and derivatives:" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 54 "In the case of a superconductor, the cost function is:
" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "f:=(x,k)->1/2*R[k]*xi[k]^
2*(dilog(exp(2*(abs(x)-j[k])/xi[k])+1)-dilog(exp(-2*j[k]/xi[k])+1))+R[
k]*xi[k]*ln(exp(2*(abs(x)-j[k])/xi[k])+1)*abs(x):" "6#>%\"fGR6$%\"xG%
\"kG7\"6$%)operatorG%&arrowG6\",&*,\"\"\"\"\"\"\"\"#!\"\"&%\"RG6#F(F1&
%#xiG6#F(\"\"#,&-%&dilogG6#,&-%$expG6#*(\"\"#F1,&-%$absG6#F'F1&%\"jG6#
F(F3F1&F86#F(F3F1\"\"\"F1F1-F=6#,&-FA6#,$*(\"\"#F1&FJ6#F(F1&F86#F(F3F3
F1\"\"\"F1F3F1F1**&F56#F(F1&F86#F(F1-%#lnG6#,&-FA6#*(\"\"#F1,&-FG6#F'F
1&FJ6#F(F3F1&F86#F(F3F1\"\"\"F1F1-FG6#F'F1F1F-F-F-" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 16 "It's derivative " }{XPPEDIT 18 0 "diff(f(x),x)" "
6#-%%diffG6$-%\"fG6#%\"xGF)" }{TEXT -1 14 " (the voltage " }{XPPEDIT 
18 0 "V(I)" "6#-%\"VG6#%\"IG" }{TEXT -1 5 ") is:" }}}{EXCHG {PARA 0 ">
 " 0 "" {XPPEDIT 19 1 "Df:=(x,k)->evalf(signum(x)*R[k]*abs(x)*(1+tanh(
(abs(x)-j[k])/xi[k]))):" "6#>%#DfGR6$%\"xG%\"kG7\"6$%)operatorG%&arrow
G6\"-%&evalfG6#**-%'signumG6#F'\"\"\"&%\"RG6#F(F5-%$absG6#F'F5,&\"\"\"
F5-%%tanhG6#*&,&-F:6#F'F5&%\"jG6#F(!\"\"F5&%#xiG6#F(FHF5F5F-F-F-" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "The second derivative " }{XPPEDIT 
18 0 "diff(f(x),`$`(x,2));" "6#-%%diffG6$-%\"fG6#%\"xG-%\"$G6$F)\"\"#
" }{TEXT -1 4 " is:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "D2f:=(
x,k)->R[k]*(1+tanh((abs(x)-j[k])/xi[k]))+R[k]*abs(x)*(1-tanh((abs(x)-j
[k])/xi[k])^2)/xi[k]:" "6#>%$D2fGR6$%\"xG%\"kG7\"6$%)operatorG%&arrowG
6\",&*&&%\"RG6#F(\"\"\",&\"\"\"F3-%%tanhG6#*&,&-%$absG6#F'F3&%\"jG6#F(
!\"\"F3&%#xiG6#F(FAF3F3F3**&F16#F(F3-F<6#F'F3,&\"\"\"F3*$-F76#*&,&-F<6
#F'F3&F?6#F(FAF3&FC6#F(FA\"\"#FAF3&FC6#F(FAF3F-F-F-" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 75 "And the most important one is the inverse of th
e voltage-current function, " }{XPPEDIT 18 0 "((`f'`)^(-1))(t);" "6#-)
%#f'G,$\"\"\"!\"\"6#%\"tG" }{TEXT -1 39 ", which involves the LambertW
 function:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "invDf:=(t,k)->e
valf(signum(t)*(1/2*abs(t)/R[k]+1/2*xi[k]*LambertW(abs(t)*exp((2*j[k]-
abs(t)/R[k])/xi[k])/(R[k]*xi[k]))) ):" "6#>%&invDfGR6$%\"tG%\"kG7\"6$%
)operatorG%&arrowG6\"-%&evalfG6#*&-%'signumG6#F'\"\"\",&**\"\"\"F5\"\"
#!\"\"-%$absG6#F'F5&%\"RG6#F(F:F5**\"\"\"F5\"\"#F:&%#xiG6#F(F5-%)Lambe
rtWG6#*(-F<6#F'F5-%$expG6#*&,&*&\"\"#F5&%\"jG6#F(F5F5*&-F<6#F'F5&F?6#F
(F:F:F5&FE6#F(F:F5*&&F?6#F(F5&FE6#F(F5F:F5F5F5F-F-F-" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 64 "The network cost function is the sum of the cos
ts over all arcs:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "F:=x->Ve
cSum(Invoke(x,f)):" "6#>%\"FGR6#%\"xG7\"6$%)operatorG%&arrowG6\"-%'Vec
SumG6#-%'InvokeG6$F'%\"fGF,F,F," }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 
"" }{TEXT 258 39 "Part III) Random values for parameters:" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 18 "The source vector:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {XPPEDIT 19 1 "ST:=vector([seq(abs(0.01*randvector(9)[i]),i=1..
n)]);" "6#>%#STG-%'vectorG6#7#-%$seqG6$-%$absG6#*&$\"\"\"!\"#\"\"\"&-%
+randvectorG6#\"\"*6#%\"iGF3/F:;\"\"\"%\"nG" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%#STG-%'vectorG6#7($\"#&)!\"#$\"#dF+$\"#aF+$\"\"\"F+$
\"#%*F+$\"#\\F+" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "Resistances, c
ritical currents and transition widths:" }}}{EXCHG {PARA 0 "> " 0 "" 
{XPPEDIT 19 1 "R:=vector([seq(abs(0.01*randvector(9)[i]),i=1..m)]);" "
6#>%\"RG-%'vectorG6#7#-%$seqG6$-%$absG6#*&$\"\"\"!\"#\"\"\"&-%+randvec
torG6#\"\"*6#%\"iGF3/F:;\"\"\"%\"mG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
#>%\"RG-%'vectorG6#7+$\"#**!\"#$\"#>F+$\"#BF+$\"#hF+$\"#<F+$\"\"%F+$\"
#iF+$\"#5F+$\"#9F+" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "j := ve
ctor([seq(abs(.1e-1*randvector(9)[k]),k = 1 .. m)]);" "6#>%\"jG-%'vect
orG6#7#-%$seqG6$-%$absG6#*&$\"\"\"!\"#\"\"\"&-%+randvectorG6#\"\"*6#%
\"kGF3/F:;\"\"\"%\"mG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"jG-%'vect
orG6#7+$\"#N!\"#$\"#VF+$\"#9F+$\"#fF+$\"#$*F+$\"#\"*F+$\"\"$F+$\"#oF+$
\"#mF+" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "xi := vector([seq(a
bs(.1e-2*randvector(9)[k]+0.025),k = 1 .. m)]);" "6#>%#xiG-%'vectorG6#
7#-%$seqG6$-%$absG6#,&*&$\"\"\"!\"$\"\"\"&-%+randvectorG6#\"\"*6#%\"kG
F4F4$\"#D!\"$F4/F;;\"\"\"%\"mG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x
iG-%'vectorG6#7+$\"#P!\"$$\"#LF+$\"#fF+$\"#WF+$\"#!*F+$\"#RF+$\"#OF+$
\"$1\"F+$\"#HF+" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 60 "Here is the ra
ndom vector of voltages for the initial guess:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {XPPEDIT 19 1 "V:=vector([seq(abs(0.01*randvector(9)[i]),i=1..n
)]);" "6#>%\"VG-%'vectorG6#7#-%$seqG6$-%$absG6#*&$\"\"\"!\"#\"\"\"&-%+
randvectorG6#\"\"*6#%\"iGF3/F:;\"\"\"%\"nG" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%\"VG-%'vectorG6#7($\"#r!\"#$\"#SF+$\"#\\F+\"\"!$\"#!*
F+$\"#?F+" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 259 54 "Part
 IV) Vectors and matrices used in the optimization" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 25 "The supply-demand vector:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {XPPEDIT 19 1 "b:=evalm(ST):" "6#>%\"bG-%&evalmG6#%#STG" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "For now, we do not assign a specif
ic value to the vector of Lagrange multipliers (voltages):" }}}{EXCHG 
{PARA 0 "> " 0 "" {XPPEDIT 19 1 "lambda:=vector(n):" "6#>%'lambdaG-%'v
ectorG6#%\"nG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 49 "The tension vect
or (potential drops across arcs):" }}}{EXCHG {PARA 0 "> " 0 "" 
{XPPEDIT 19 1 "t:=multiply(transpose(A),lambda):" "6#>%\"tG-%)multiply
G6$-%*transposeG6#%\"AG%'lambdaG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
72 "The flow vector (currents) for the given potentials can now be fou
nd as " }{XPPEDIT 18 0 "x[`*`]=(`f'`^(-1))(t)" "6#/&%\"xG6#%\"*G-)%#f'
G,$\"\"\"!\"\"6#%\"tG" }{TEXT -1 2 ": " }}}{EXCHG {PARA 0 "> " 0 "" 
{XPPEDIT 19 1 "x[`*`]:=Invoke(t,invDf):" "6#>&%\"xG6#%\"*G-%'InvokeG6$
%\"tG%&invDfG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "And the gradient
 is:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "G:=evalm(multiply(A,x
[`*`])-b):" "6#>%\"GG-%&evalmG6#,&-%)multiplyG6$%\"AG&%\"xG6#%\"*G\"\"
\"%\"bG!\"\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 98 "The diagonal of t
he conjugate Hessian is the inverse of the diagonal of the cost-functi
on Hessian:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "DH[`*`]:=Diago
nal(Invoke(x[`*`],(x,i)->1/D2f(x,i))):" "6#>&%#DHG6#%\"*G-%)DiagonalG6
#-%'InvokeG6$&%\"xG6#F'R6$F/%\"iG7\"6$%)operatorG%&arrowG6\"*&\"\"\"\"
\"\"-%$D2fG6$F/F3!\"\"F8F8F8" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "A
nd the full Hessian is:" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "H:
=multiply(A,DH[`*`],transpose(A)):" "6#>%\"HG-%)multiplyG6%%\"AG&%#DHG
6#%\"*G-%*transposeG6#F(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }
{TEXT 260 41 "Part IV) Taking a step toward the minimum" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT 261 76 "This piece of code would be repeated unt
il convergence to the global minimum" }{TEXT -1 1 "." }}{PARA 0 "" 0 "
" {TEXT -1 88 "We start by assigning a random vector of voltages to th
e vector of Lagrange multipliers:" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "
lambda:=evalm(V):" "6#>%'lambdaG-%&evalmG6#%\"VG" }}{PARA 0 "" 0 "" 
{TEXT -1 58 "The direction vector is the solution to the Newton system
 " }{XPPEDIT 18 0 "H*d=-G" "6#/*&%\"HG\"\"\"%\"dGF&,$%\"GG!\"\"" }
{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "d:=linsolve(map(eval
,H),evalm(-G)):" "6#>%\"dG-%)linsolveG6$-%$mapG6$%%evalG%\"HG-%&evalmG
6#,$%\"GG!\"\"" }}{PARA 0 "" 0 "" {TEXT -1 45 "Now we do a line search
 along this direction " }{XPPEDIT 18 0 "lambda=lambda+alpha*d" "6#/%'l
ambdaG,&F$\"\"\"*&%&alphaGF&%\"dGF&F&" }{TEXT -1 1 ":" }}{PARA 0 "> " 
0 "" {XPPEDIT 19 1 "lambda := evalm(V+alpha*d):" "6#>%'lambdaG-%&evalm
G6#,&%\"VG\"\"\"*&%&alphaGF*%\"dGF*F*" }}{PARA 0 "" 0 "" {TEXT -1 41 "
The objective function is the Lagrangian " }{XPPEDIT 18 0 "h(alpha)=L[
lambda+alpha*d]" "6#/-%\"hG6#%&alphaG&%\"LG6#,&%'lambdaG\"\"\"*&F'F-%
\"dGF-F-" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "L[lambda]=x[`*`]^T*t-
F(x[`*`])-lambda^T*d" "6#/&%\"LG6#%'lambdaG,(*&)&%\"xG6#%\"*G%\"TG\"\"
\"%\"tGF0F0-%\"FG6#&F,6#F.!\"\"*&)F'F/F0%\"dGF0F7" }{TEXT -1 1 ":" }}
{PARA 0 "> " 0 "" {XPPEDIT 19 1 "L[lambda]:=unapply(dotprod(map(eval,x
[`*`]),map(eval,t),'orthogonal')-F(map(eval,x[`*`]))-dotprod(lambda,b,
'orthogonal'),alpha):" "6#>&%\"LG6#%'lambdaG-%(unapplyG6$,(-%(dotprodG
6%-%$mapG6$%%evalG&%\"xG6#%\"*G-F06$F2%\"tG.%+orthogonalG\"\"\"-%\"FG6
#-F06$F2&F46#F6!\"\"-F-6%F'%\"bG.F;FD%&alphaG" }}{PARA 0 "" 0 "" 
{TEXT -1 48 "This function should be concave and continious: " }}
{PARA 0 "> " 0 "" {XPPEDIT 19 1 "plot(L[lambda](alpha),alpha=0..1,nump
oints=10,adaptive=false,labels=[\"alpha\",\"L[lambda](alpha)\"]);" "6#
-%%plotG6'-&%\"LG6#%'lambdaG6#%&alphaG/F,;\"\"!\"\"\"/%*numpointsG\"#5
/%)adaptiveG%&falseG/%'labelsG7$Q&alpha6\"Q1L[lambda](alpha)F;" }}
{PARA 0 "" 0 "" {TEXT -1 75 "To find the minimum (i.e. perform the lin
e search), we find the derivative " }{XPPEDIT 18 0 "h(alpha)=diff(L[la
mbda](alpha),alpha)" "6#/-%\"hG6#%&alphaG-%%diffG6$-&%\"LG6#%'lambdaG6
#F'F'" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "h(alpha)=d^T*G" "6#/-%\"
hG6#%&alphaG*&)%\"dG%\"TG\"\"\"%\"GGF," }{TEXT -1 21 ", and set it to \+
zero:" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "h:=unapply(dotprod(d,map(eva
l,G),'orthogonal'),alpha):" "6#>%\"hG-%(unapplyG6$-%(dotprodG6%%\"dG-%
$mapG6$%%evalG%\"GG.%+orthogonalG%&alphaG" }}{PARA 0 "" 0 "" {TEXT -1 
166 "The plot below shows that h(alpha) is not continious. But this is
 an optical illusion, it comes from the fact that there is a very shar
p non-linear transition in the " }{XPPEDIT 18 0 "V(I)" "6#-%\"VG6#%\"I
G" }{TEXT -1 40 " curve. The derivative of this function " }{XPPEDIT 
18 0 "diff(h(alpha),alpha)=d^T*H*d" "6#/-%%diffG6$-%\"hG6#%&alphaGF**(
)%\"dG%\"TG\"\"\"%\"HGF/F-F/" }{TEXT -1 88 " is this practically undef
ined (infinite) at certain points. The theory guarantees that " }
{XPPEDIT 18 0 "h(alpha)" "6#-%\"hG6#%&alphaG" }{TEXT -1 202 " is piece
wise convex, continious, and non-decreasing. So Newton's method for th
e line search minimization, if used, must be used with great care and \+
in combination with secant-like or bisection methods!" }}{PARA 0 "> " 
0 "" {XPPEDIT 19 1 "dh:=unapply(dotprod(d,multiply(map(eval,H),d),'ort
hogonal'),alpha):" "6#>%#dhG-%(unapplyG6$-%(dotprodG6%%\"dG-%)multiply
G6$-%$mapG6$%%evalG%\"HGF+.%+orthogonalG%&alphaG" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "plot(h(alpha),alpha =
 0 .. 1.5,labels=[\"alpha\",\"h(alpha)\"]);" "6#-%%plotG6%-%\"hG6#%&al
phaG/F);\"\"!$\"#:!\"\"/%'labelsG7$Q&alpha6\"Q)h(alpha)F4" }}{PARA 13 
"" 1 "" {GLPLOT2D 419 206 206 {PLOTDATA 2 "6%-%'CURVESG6$7,7$\"\"!$\"+
F@pZ>!\"*7$$\"+*)[^i6!#5$\"+26ZT7F+7$$\"+yG,u@F/$\"+*>!p*Q)F/7$$\"+b\"
[:J$F/$\"+O7$)*o'F/7$$\"+bPicWF/$\"+i*>?k'F/7$$\"+6qD'f&F/$\"+ZW8E\")F
/7$$\"+AH%Gl'F/$\"+9gk\"4\"F+7$$\"+Lk(ou(F/$\"+mX#yc\"F+7$$\"+buKy))F/
$\"+j.M5AF+7$$\"\"\"F($\"+pg1mIF+-%'COLOURG6&%$RGBG$\"#5!\"\"F(F(-%+AX
ESLABELSG6$Q&alpha6\"Q1L[lambda](alpha)F_o-%%VIEWG6$;F(FV%(DEFAULTG" 
1 2 0 1 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 }}}{PARA 13 "" 1 "
" {GLPLOT2D 416 214 214 {PLOTDATA 2 "6%-%'CURVESG6$7^q7$\"\"!$!+_Jy\\q
!\"*7$$\"+DJdpK!#6$!+%H1u\"oF+7$$\"+t@*>p%F/$!+1y#zq'F+7$$\"+>7T9hF/$!
+<VYvlF+7$$\"+mfK9lF/$!+AwPClF+7$$\"+72C9pF/$!+;]yYkF+7$$\"+&3)>9rF/$!
+:+_ajF+7$$\"+ea:9tF/$!+JqJ8_F+7$$\"+JG69vF/$!+^q@Y^F+7$$\"+/-29xF/$!+
p.H1^F+7$$\"+'p**Q^)F/$!+`#>\\+&F+7$$\"+)=HPJ*F/$!+=<HI\\F+7$$\"+v\"*R
#4\"!#5$!+TKz,[F+7$$\"+JaU`7Fao$!+5no#o%F+7$$\"+%GZRd\"Fao$!+**H^VWF+7
$$\"+y'HDs\"Fao$!+JPH?VF+7$$\"+s?6r=Fao$!+e\")HoTF+7$$\"+IJM!*=Fao$!+D
dLUTF+7$$\"+)=u&4>Fao$!+n4u6TF+7$$\"+<(*=>>Fao$!+GM3$4%F+7$$\"+Y_!)G>F
ao$!+D?dpSF+7$$\"+v2UQ>Fao$!+F*QE.%F+7$$\"+/j.[>Fao$!+d!3Kh#F+7$$\"+it
En>Fao$!+\")HEBDF+7$$\"+?%)\\')>Fao$!+n!QV[#F+7$$\"+y%Hd+#Fao$!+bJO`CF
+7$$\"+N0'\\-#Fao$!+XBeDCF+7$$\"+nZ)=5#Fao$!+xDgCBF+7$$\"+(**3)y@Fao$!
+LiGEAF+7$$\"+)H>zL#Fao$!+G,@>?F+7$$\"+(fHq\\#Fao$!+O0$3\"=F+7$$\"+f'H
U\"GFao$!+A3v*R\"F+7$$\"+8*309$Fao$!+*eDf\")*Fao7$$\"+ce*yU$Fao$!+K,Za
hFao7$$\"+)[D9v$Fao$!+>NjY?Fao7$$\"+jNGwSFao$\"+tD<c?Fao7$$\"+8XM*Q%Fa
o$\"+++7sdFao7$$\"+ZQjtYFao$\"+^&eJ-*Fao7$$\"+j8o6]Fao$\"+[ED'G\"F+7$$
\"+]>0)H&Fao$\"+(yE7h\"F+7$$\"+A!p6j&Fao$\"+_K=!*>F+7$$\"+vS.EfFao$\"+
SgoFBF+7$$\"+C4z(3'Fao$\"+%z(R9DF+7$$\"+sxa\\iFao$\"+#p/Kq#F+7$$\"+eJc
EjFao$\"+Vu^%z#F+7$$\"+W&yNS'Fao$\"+rGW))GF+7$$\"+\"RKGU'Fao$\"+qt=8HF
+7$$\"+Pi3UkFao$\"+kJ()RHF+7$$\"+hJr^kFao$\"+;1obHF+7$$\"+%3S8Y'Fao$\"
+&[g)*\\$F+7$$\"+2q'4Z'Fao$\"+kzBFNF+7$$\"+IRf!['Fao$\"+5SXUNF+7$$\"+B
;5>lFao$\"+*e@Kf$F+7$$\"+;$4wb'Fao$\"+1y4SOF+7$$\"+oUK=nFao$\"+f>kDQF+
7$$\"+>#R!zoFao$\"+TcM0SF+7$$\"+4A@urFao$\"+iDKOVF+7$$\"+chf#\\(Fao$\"
+;8r-ZF+7$$\"+(f2L#yFao$\"+H5,*3&F+7$$\"+yG>6\")Fao$\"+5Z,FaF+7$$\"+po
6A%)Fao$\"+%3(*[z&F+7$$\"+BVs#e)Fao$\"+Fm*y)fF+7$$\"+v<LV()Fao$\"+#y:o
='F+7$$\"+1L*=#))Fao$\"+Tv%*)G'F+7$$\"+N[X+*)Fao$\"+[#4jR'F+7$$\"+,ctR
*)Fao$\"+'y'R`kF+7$$\"+lj,z*)Fao$\"+NZ[>lF+7$$\"+dl$)))*)Fao$\"+\"RiWa
'F+7$$\"+[nl)**)Fao$\"+^r,xqF+7$$\"+RpZ3!*Fao$\"+1=a/rF+7$$\"+IrH=!*Fa
o$\"+'=PP7(F+7$$\"+7v$z.*Fao$\"+=7YcrF+7$$\"+%*ydd!*Fao$\"+0GQ'=(F+7$$
\"+LDg4#*Fao$\"+v\"\\wR(F+7$$\"+srih$*Fao$\"+MNZ+wF+7$$\"+c?A*p*Fao$\"
+0,gW!)F+7$$\"+3mD+5F+$\"+xLQT%)F+7$$\"+c]kK5F+$\"+Pjhk))F+7$$\"+0Q*>1
\"F+$\"+(G_([#*F+7$$\"+R(zS4\"F+$\"+;X!3n*F+7$$\"+-,FC6F+$\"+E;i15!\")
7$$\"+Jx#e:\"F+$\"+;YuZ5Fhdl7$$\"+\"3\"o'=\"F+$\"+'*zf)3\"Fhdl7$$\"+!o
\")*=7F+$\"+]KFJ6Fhdl7$$\"+&*44]7F+$\"+\"G2B<\"Fhdl7$$\"+jZ!>G\"F+$\"+
\"ogU@\"Fhdl7$$\"+(4bMJ\"F+$\"+k:!fD\"Fhdl7$$\"+ylWU8F+$\"+68C%H\"Fhdl
7$$\"+'3ucP\"F+$\"+W+PQ8Fhdl7$$\"+lJR09F+$\"+,)Q#y8Fhdl7$$\"+MlB@9F+$
\"+1=#**R\"Fhdl7$$\"+-*zqV\"F+$\"+nTFB9Fhdl7$$\"+l>mW9F+$\"+IIl3:Fhdl7
$$\"+GSC_9F+$\"+ly#)\\:Fhdl7$$\"+\"4E)f9F+$\"+(e,Ac\"Fhdl7$$\"+`\"3uY
\"F+$\"+ZM.t:Fhdl7$$\"+xSq$[\"F+$\"+OeV&f\"Fhdl7$$\"#:!\"\"$\"+&H2uh\"
Fhdl-%'COLOURG6&%$RGBG$\"#5F\\jlF(F(-%+AXESLABELSG6$Q&alpha6\"Q)h(alph
a)Fijl-%%VIEWG6$;F(Fjil%(DEFAULTG" 1 2 0 1 0 2 9 1 4 2 1.000000 
45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "Now we f
ind the zero of " }{XPPEDIT 18 0 "h(alpha)" "6#-%\"hG6#%&alphaG" }
{TEXT -1 34 "--Beware, fsolve often fails here:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {XPPEDIT 19 1 "alpha:=fsolve(h(alpha)=0,alpha=0..2);" "6#>%&alp
haG-%'fsolveG6$/-%\"hG6#F$\"\"!/F$;F,\"\"#" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%&alphaG$\"+K]$H\"R!#5" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 43 "And the new estimate for the potentials is:" }}}{EXCHG 
{PARA 0 "> " 0 "" {XPPEDIT 19 1 "V:=map(eval,lambda);" "6#>%\"VG-%$map
G6$%%evalG%'lambdaG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'vectorG6#7($
\"+MgRHp!#5$\"+>IT15!\"*$\"+(Gq=1'F)$\"+$3fo(>!#7$\"+y3#Hd'F)$\"+@Y<=E
F)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "This of course has not conv
erged yet." }}}}{MARK "1 0 2" 2 }{VIEWOPTS 1 1 0 1 1 1803 }