content

• Preface
• Orthogonal Transformation and Least Squares
• Cholesky Decomposition
• The Square Root Free Cholesky Algorithm
• The Householder Transformation
• Modified Gram–Schmidt orthogonalization
• Solving nolinear Equation
• Newton-Raphson Method
• mulit-root iterative improvement
• Muhammad method
• numerical Optimization
• Successive parabolic interpolation
• Newton’s Method
• Quasi-Newton Methods(DFP)
• Quasi-Newton Methods(BFGS)
• Interpolation
• Lagrange interpolation
• Newton method
• Hermit interpolation
• data fitting
• Polynomial fitting ephemeris
• Chebyshev fitting ephemeris
• ODE
• Symplectic Geometric Algorithms for Hamiltonian Systems
• Precise integration method
• Appendix：Matrix and Linear algebra class design
• C++ matrix class design
• C# matrix class design
• golang's slice for linear algebra

数值优化

一维搜索之连续抛物线插值

在黄金分割法中，除了比较函数值大小之外，函数值本身并未在最优求解中用到。 这说明，我们并未充分利用现有的信息来加速迭代收敛。这一节介绍连续抛物线插值法 求解一维搜索问题，这是一种较少浪费函数值的方法。
下面阐述连续抛物线插值的基本思想，先在最小值附件选3个点$r,s,t$。计算这3个点的目标函数值$f$，且可以 给出过这3个点的抛物线 $$P\left(x\right)=f\left(r\right)+d_1\left(x-r\right)+d_3\left(x-r\right)\left(x-s\right)$$ 其中 $$d_1=\frac{f\left(s\right)-f\left(r\right)}{s-r}$$ $$d_2=\frac{f\left(t\right)-f\left(s\right)}{t-s}$$ $$d_3=\frac{d_2-d_1}{t-r}$$ 令$P\left(x\right)$可以给出抛物线的极小值公式 $$x=\frac{r+s}{2}-\frac{\left(f\left(s\right)-f\left(r\right)\right)\left(t-r\right)\left(t-s\right)}{2\left[\left(s-t\right)\left(f\left(t\right)-f\left(s\right)\right)-\left(f\left(s\right)-f\left(r\right)\right)\left(t-s\right)\right]}$$ 此点即为极小值新的近似。有些教材给的算法到此为止，这样的话只能给出较粗的极小值。 在实际应用时，可以重复以上过程，用新的x代替r，s，t中最远的或次最远的点， 反复迭代，即可以给出需求精度下的极小值。
下面给出连续抛物线求一维搜索算法：

• [Input] 最大允许迭代次数、误差容限，初始猜想极小点r，s，t 。
• [Output] 极值点的自变量与函数值。
• [STEP1] 对$i=1,2,3,\cdots$，作如下循环 $$x=\frac{r+s}{2}-\frac{\left(f\left(s\right)-f\left(r\right)\right)\left(t-r\right)\left(t-s\right)}{2\left[\left(s-t\right)\left(f\left(t\right)-f\left(s\right)\right)-\left(f\left(s\right)-f\left(r\right)\right)\left(t-s\right)\right]}$$ $$t=s$$ $$s=r$$ $$r=x$$ 如果$i$超过与设定的最大允许迭代次数，停止循环，并给出告警信息。 如果已经满足预设的精度需求，则停止循环。
• [STEP2] 输出计算结果，并停止该算法。

一般而言，连续抛物线插值法要快于黄金分割方法，因为它更有效的利用了函数值的信息。 用连续抛物线插值法求解函数 $$f\left(x\right)=x^6-11x^3+17x^2-7x+1$$ 在区间$\left[0,1\right]$上的极小值。 go语言版本代码如下。

package main
import (
	"fmt"
	."math"
)
func main() {
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-7 18:25
  go build -gcflags "-N -l"  
     
------------------------------------------------------*/
	r :=0.0
	s :=0.7
	t :=1.0
	_ , _ = paraInterp(r,s,t,myfx)

}

func myfx(x float64) float64 {
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-7 18:26
  go build -gcflags "-N -l"  
------------------------------------------------------*/
	fx := Pow(x,6)-11.0*Pow(x,3)+17.0*x*x - 7.0*x +1.0
	return fx
}

func paraInterp(r,s,t float64,funcX func(float64) float64) (x, fx float64) {
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-7 18:28
  go build -gcflags "-N -l"  
     连续抛物型法一维搜索
  -----------------------------------------------------  
  Input  Parameters :
     r,s,t  ------初值
     funcX -----  目标函数
  Output Parameters :
     x---极小点
     fx--极小点函数值
------------------------------------------------------*/
	var imax int = 20
	var tol float64 = 1e-8
	
	var fs,ft,fr float64
	var tmp1,tmp2 float64
	
	fs = funcX(s)
	fr = funcX(r)
	ft = funcX(t)

	for i := 0; i < imax; i++ {		
		tmp1 = (fs -fr) * (t-r) * (t-s)
		tmp2 = 2.0*((s-r)*(ft-fs) -(fs-fr)*(t-s))
		
		x = (r+s)/2.0 - tmp1/tmp2
		
		if Abs(fs-fr) < tol {
		    break
		}
		
	    t = s
		s = r
		r = x
		
		ft = fs
		fs = fr
		fr = funcX(r)
		// *+++single function evaluation +++*   
		fmt.Printf("i= %4d  x= %12.7f  f(x)=%12.7f \n", i, x, ft)
	}
	return x,fx
}