// ODEMilneSimpson
/*
------------------------------------------------------
作者   : Black Ghost
日期   : 2018-12-26
版本   : 0.0.0
------------------------------------------------------
    Milne-Simpson预估校正方法
理论：
    预估：
                         4h
    p_(k+1) = y_(k-3) + ---(2f_(k-2)-f_(k-1)+2fk)
                         3
    校正：
                         h
    y_(k+1) = y_(k-1) + ---(f_(k-1)+4fk+f_(k+1))
                         3

    步长 h < 0.45/|fy(x,y)|

    四阶精度

    参考：John H. Mathews and Kurtis D. Fink. Numerical
         methods using MATLAB, 4th ed. Pearson
         Education, 2004. ss 9.6.4
------------------------------------------------------
输入   :
    fun     被积分函数
    x0      初值,2x4
    h       步长
    n       积分步数
输出   :
    sol     解矩阵
    err     解出标志：false-未解出或达到步数上限；
                     true-全部解出
------------------------------------------------------
*/

package goNum

// ODEMilneSimpson Milne-Simpson预估校正方法
func ODEMilneSimpson(fun func(float64, float64) float64, x0 Matrix, h float64, n int) (Matrix, bool) {
	/*
		Milne-Simpson预估校正方法
		输入   :
		    fun     被积分函数
		    x0      初值,2x4
		    h       步长
		    n       积分步数
		输出   :
		    sol     解矩阵
		    err     解出标志：false-未解出或达到步数上限；
		                     true-全部解出
	*/
	//判断n
	if n < 0 {
		panic("Error in goNum.ODEMilneSimpson: n is not a positive value")
	}
	//判断初值
	if (x0.Rows != 2) || (x0.Columns < 4) {
		panic("Error in goNum.ODEMilneSimpson: Initial values error")
	}

	sol := ZeroMatrix(2, n+1)
	p := ZeroMatrix(n+1, 1)
	var err bool = false

	//初值
	for i := 0; i < 4; i++ {
		sol.SetMatrix(0, i, x0.GetFromMatrix(0, i))
		sol.SetMatrix(1, i, x0.GetFromMatrix(1, i))
	}

	//计算
	for i := 4; i < n+1; i++ {
		sol.SetMatrix(0, i, sol.GetFromMatrix(0, i-1)+h) //xi
		//pi
		temp0 := fun(sol.GetFromMatrix(0, i-3), sol.GetFromMatrix(1, i-3)) //f_(i-3)
		temp1 := fun(sol.GetFromMatrix(0, i-2), sol.GetFromMatrix(1, i-2)) //f_(i-2)
		temp2 := fun(sol.GetFromMatrix(0, i-1), sol.GetFromMatrix(1, i-1)) //f_(i-1)
		soltemp := 2.0 * temp0
		soltemp += -1.0 * temp1
		soltemp += 2.0 * temp2
		p.SetMatrix(i, 0, sol.GetFromMatrix(1, i-4)+4.0*h*soltemp/3.0)
		//yi
		soltemp = temp1
		soltemp += 4.0 * temp2
		soltemp += fun(sol.GetFromMatrix(0, i), p.GetFromMatrix(i, 0)) //fi
		sol.SetMatrix(1, i, sol.GetFromMatrix(1, i-2)+h*soltemp/3.0)
	}

	err = true
	return sol, err
}