Milne格式MATLAB,ODEMilneSimpson.go

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

}

一键复制

编辑

Web IDE

原始数据

按行查看

历史