mirror of
https://gitee.com/chfenger/goNum.git
synced 2025-12-06 08:48:55 +08:00
96fccdb6ef
update comments in all of files in convenient to generate package information using `godoc` command or give tips in LiteIDE editor when the cursor keeps resting on a function or method.
195 lines
5.9 KiB
Go
195 lines
5.9 KiB
Go
// InterpSpline22
|
||
/*
|
||
------------------------------------------------------
|
||
作者 : Black Ghost
|
||
日期 : 2018-12-9
|
||
版本 : 0.0.0
|
||
------------------------------------------------------
|
||
用节点处的二阶导数表示的三次样条插值函数,
|
||
二阶导数边界条件
|
||
n+1个点, n个区间
|
||
理论：
|
||
区间[x(i-1), xi]上的三次样条函数表达为：
|
||
(xi-x)^3
|
||
Si(x) = ----------M(i-1) +
|
||
6*hi
|
||
(x-x(i-1))^3
|
||
--------------Mi +
|
||
6*hi
|
||
M(i-1) xi-x
|
||
(y(i-1) - --------hi^2)-------
|
||
6 hi
|
||
Mi x-x(i-1)
|
||
(yi - ----hi^2)----------
|
||
6 hi
|
||
|
||
令 Mi = hi/(hi+h(i+1))
|
||
lambdai = 1-Mi = h(i+1)/(hi+h(i+1))
|
||
6 y(i+1)-yi yi-y(i-1)
|
||
fi = -----------(---------- - -----------)
|
||
hi+h(i+1) h(i+1) hi
|
||
(i = 1,...,n-1)
|
||
|
||
则mi可由n-1阶线性方程组求得（利用LEs_Chasing）：
|
||
|2 l1 || M1 | | f1-M1*M0 |
|
||
| M2 2 l2 || M2 | = | f2 |
|
||
| ........ || ... | | ... |
|
||
| M(n-2) 2 l(n-2)||M(n-2)| | f(n-2) |
|
||
| M(n-1) 2 ||M(n-1)| |f(n-1)-l(n-1)Mn|
|
||
|
||
参考 李信真, 车刚明, 欧阳洁, 等. 计算方法. 西北工业大学
|
||
出版社, 2000, pp 124-127.
|
||
------------------------------------------------------
|
||
输入 :
|
||
A 数据点矩阵，(n+1)x3，第一列xi；第二列yi；
|
||
第三列y''i，且y''i只需给出y''0和y''n
|
||
输出 :
|
||
B 插值方程系数结果矩阵，从前到后对应从0到3阶，4xn
|
||
err 解出标志：false-未解出或达到步数上限；
|
||
true-全部解出
|
||
------------------------------------------------------
|
||
*/
|
||
|
||
package goNum
|
||
|
||
// InterpSpline22 用节点处的二阶导数表示的三次样条插值函数, 二阶导数边界条件
|
||
func InterpSpline22(A Matrix) (Matrix, bool) {
|
||
/*
|
||
用节点处的二阶导数表示的三次样条插值函数, 二阶导数边界条件
|
||
输入 :
|
||
A 数据点矩阵，(n+1)x3，第一列xi；第二列yi；
|
||
第三列y'i，且y'i只需给出y'0和y'n
|
||
输出 :
|
||
B 插值方程系数结果矩阵，从前到后对应从0到3阶，4xn
|
||
err 解出标志：false-未解出或达到步数上限；
|
||
true-全部解出
|
||
*/
|
||
var err bool = false
|
||
n := A.Rows - 1
|
||
sol := ZeroMatrix(4, n)
|
||
BA := ZeroMatrix(n-1, n-1) //对角占优的三对角矩阵
|
||
BB := ZeroMatrix(n-1, 1) //解向量
|
||
BC := ZeroMatrix(n-1, 1) //值向量
|
||
|
||
//1解插值函数的一阶导数mi
|
||
//1.0.1第一行
|
||
if true { //限制变量使用范围
|
||
h1 := A.GetFromMatrix(1, 0) - A.GetFromMatrix(0, 0)
|
||
h2 := A.GetFromMatrix(2, 0) - A.GetFromMatrix(1, 0)
|
||
y0 := A.GetFromMatrix(0, 1)
|
||
y1 := A.GetFromMatrix(1, 1)
|
||
y2 := A.GetFromMatrix(2, 1)
|
||
M1 := h1 / (h1 + h2)
|
||
l1 := 1.0 - M1
|
||
f1 := 6.0 * ((y2-y1)/h2 - (y1-y0)/h1) / (h1 + h2)
|
||
BA.SetMatrix(0, 0, 2.0)
|
||
BA.SetMatrix(0, 1, l1)
|
||
BC.Data[0] = f1 - M1*A.GetFromMatrix(0, 2)
|
||
}
|
||
//1.0.2其它行
|
||
for i := 2; i < n-1; i++ {
|
||
yi_1 := A.GetFromMatrix(i-1, 0)
|
||
yi := A.GetFromMatrix(i, 0)
|
||
yi1 := A.GetFromMatrix(i+1, 0)
|
||
|
||
hi := A.GetFromMatrix(i, 0) - A.GetFromMatrix(i-1, 0)
|
||
hi1 := A.GetFromMatrix(i+1, 0) - A.GetFromMatrix(i, 0)
|
||
|
||
Mi := hi / (hi + hi1)
|
||
li := 1.0 - Mi
|
||
fi := 6.0 * ((yi1-yi)/hi1 - (yi-yi_1)/hi) / (hi + hi1)
|
||
//赋予BA
|
||
BA.SetMatrix(i-1, i-2, Mi)
|
||
BA.SetMatrix(i-1, i-1, 2.0)
|
||
BA.SetMatrix(i-1, i, li)
|
||
BC.Data[i-1] = fi
|
||
}
|
||
//1.0.3最后一行
|
||
if true { //i=n-1
|
||
hn_1 := A.GetFromMatrix(n-1, 0) - A.GetFromMatrix(n-2, 0)
|
||
hn := A.GetFromMatrix(n, 0) - A.GetFromMatrix(n-1, 0)
|
||
yn_2 := A.GetFromMatrix(n-2, 1)
|
||
yn_1 := A.GetFromMatrix(n-1, 1)
|
||
yn := A.GetFromMatrix(n, 1)
|
||
|
||
Mn_1 := hn_1 / (hn_1 + hn)
|
||
ln_1 := 1.0 - Mn_1
|
||
fn_1 := 6.0 * ((yn-yn_1)/hn - (yn_1-yn_2)/hn_1) / (hn_1 + hn)
|
||
|
||
BA.SetMatrix(n-2, n-3, Mn_1)
|
||
BA.SetMatrix(n-2, n-2, 2.0)
|
||
BC.Data[n-2] = fn_1 - ln_1*A.GetFromMatrix(n, 2)
|
||
}
|
||
//1.1求解
|
||
soltemp, errtemp := LEs_Chasing(BA, BC)
|
||
if errtemp != true {
|
||
panic("Error in goNum.InterpSpline11: Solve Error with goNum.LEs_Chasing")
|
||
}
|
||
for i := 0; i < n-1; i++ {
|
||
BB.Data[i] = soltemp.Data[i]
|
||
}
|
||
|
||
//2求解Si(x)
|
||
S0 := ZeroMatrix(4, 1)
|
||
S1 := ZeroMatrix(4, 1)
|
||
S2 := ZeroMatrix(4, 1)
|
||
S3 := ZeroMatrix(4, 1)
|
||
for i := 1; i < n+1; i++ {
|
||
xi_1 := A.GetFromMatrix(i-1, 0)
|
||
xi := A.GetFromMatrix(i, 0)
|
||
yi_1 := A.GetFromMatrix(i-1, 1)
|
||
yi := A.GetFromMatrix(i, 1)
|
||
Mi_1 := 0.0
|
||
Mi := 0.0
|
||
if i == 1 {
|
||
Mi_1 = A.GetFromMatrix(0, 2)
|
||
Mi = BB.Data[i-1]
|
||
} else if i == n {
|
||
Mi_1 = BB.Data[i-2]
|
||
Mi = A.GetFromMatrix(n, 2)
|
||
} else {
|
||
Mi_1 = BB.Data[i-2]
|
||
Mi = BB.Data[i-1]
|
||
}
|
||
hi := xi - xi_1
|
||
temp0 := ZeroMatrix(4, 1)
|
||
//2.1 S0
|
||
temp0.Data[3] = -1.0
|
||
temp0.Data[2] = 3.0 * xi
|
||
temp0.Data[1] = -3.0 * xi * xi
|
||
temp0.Data[0] = xi * xi * xi
|
||
for j := 0; j < 4; j++ {
|
||
S0.Data[j] = temp0.Data[j] * Mi_1 / (6.0 * hi)
|
||
}
|
||
//2.1 S1
|
||
temp0.Data[3] = 1.0
|
||
temp0.Data[2] = -3.0 * xi_1
|
||
temp0.Data[1] = 3.0 * xi_1 * xi_1
|
||
temp0.Data[0] = -1.0 * xi_1 * xi_1 * xi_1
|
||
for j := 0; j < 4; j++ {
|
||
S0.Data[j] = temp0.Data[j] * Mi / (6.0 * hi)
|
||
}
|
||
//2.2 S2
|
||
temp0 = ZeroMatrix(4, 1)
|
||
temp0.Data[1] = -1.0
|
||
temp0.Data[0] = xi
|
||
for j := 0; j < 4; j++ {
|
||
S2.Data[j] = temp0.Data[j] * (yi_1 - Mi_1*hi*hi/6.0) / hi
|
||
}
|
||
//2.3 S3
|
||
temp0 = ZeroMatrix(4, 1)
|
||
temp0.Data[1] = 1.0
|
||
temp0.Data[0] = -1.0 * xi_1
|
||
for j := 0; j < 4; j++ {
|
||
S3.Data[j] = temp0.Data[j] * (yi - Mi*hi*hi/6.0) / hi
|
||
}
|
||
//2.4 Si(x)
|
||
for j := 0; j < 4; j++ {
|
||
sol.SetMatrix(j, i-1, S0.Data[j]+S1.Data[j]+S2.Data[j]+S3.Data[j])
|
||
}
|
||
}
|
||
|
||
err = true
|
||
return sol, err
|
||
}
|