代码从以下书本中的相关章节改写而来,为简化问题,这里只用来拟合三次多项式函数。
M3 的代码是列主元消去法解线性方程:
import scala.math._
object M3 {
def eliminate(ac: Array[Array[Double]], bc: Array[Double], k: Int): Unit = {
val n = bc.length - 1
var m: Double = 0.1
for (i <- k + 1 to n) {
m = ac(i)(k) / ac(k)(k)
for (j <- k + 1 to n) {
ac(i)(j) = ac(i)(j) - m * ac(k)(j)
}
bc(i) = bc(i) - m * bc(k)
//println(bc(i))
}
}
def backSubstitute(ac: Array[Array[Double]], bc: Array[Double]): Unit = {
val n = bc.length - 1
bc(n) = bc(n) / ac(n)(n)
for (i <- n - 1 to 0 by -1) {
var s = 0.0
for (j <- i+1 to n) {
s = s + ac(i)(j) * bc(j)
}
bc(i) = (bc(i) - s) / ac(i)(i)
}
}
def selectColumnPivot(ac: Array[Array[Double]], k: Int): (Int, Double) = {
var ik: Int = 0
var acMax: Double = 0.0
val n = ac.length - 1
for (i <- k to n) {
if (abs(ac(i)(k)) > acMax) {
acMax = abs(ac(i)(k))
ik = i
}
}
(ik, acMax)
}
def columnPivotGauss(ac: Array[Array[Double]], bc: Array[Double]): Unit = {
val n = bc.length - 1
var ik = 0
var acMax = 0.00
for (k <- 0 to n-1){
var (ik, acMax) = selectColumnPivot(ac, k)
if (ik != k) {
for (j <- k to n) {
var sa = (ac(k)(j), ac(ik)(j))
ac(k)(j) = sa._2
ac(ik)(j) = sa._1
}
var sb = (bc(k), bc(ik))
bc(k) = sb._2
bc(ik) = sb._1
}
eliminate(ac, bc, k)
}
backSubstitute(ac, bc)
}
}
M2代码用于把线性拟合问题转化为解线性方程组的问题:
import scala.math._
object M2 {
def fxi(x: Double, m: Int): Double = {
if (m == 0) 1.0 else pow(x, m)
}
def vectorATimesB(ac: Array[Double], bc: Array[Double], wc: Array[Double]): Double = {
val n = ac.length
(0 until n).map(i => ac(i) * bc(i) * wc(i)).sum
}
def getVI(xc: Array[Double], m: Int): Array[Double] = {
xc.map(i => fxi(i, m)).toArray
}
def formAandB(xc: Array[Double], yc: Array[Double], wc: Array[Double]):
(Array[Array[Double]], Array[Double]) = {
val xl = xc.length - 1
val n = 4
var ma: Array[Array[Double]] = Array.ofDim[Double](4, 4)
var mb: Array[Double] = new Array(4)
for (cm <- 0 to 3) {
var vi: Array[Double] = getVI(xc, cm)
ma(cm)(cm) = vectorATimesB(vi, vi, wc)
mb(cm) = vectorATimesB(vi, yc, wc)
for (j <- cm+1 to 3) {
var vj: Array[Double] = getVI(xc, j)
ma(cm)(j) = vectorATimesB(vi, vj, wc)
ma(j)(cm) = ma(cm)(j)
}
}
(ma, mb)
}
def sqrFit(xc: Array[Double], yc: Array[Double], wc: Array[Double]): Array[Double] = {
val AandB = formAandB(xc, yc, wc)
val ma = AandB._1
var mb = AandB._2
M3.columnPivotGauss(ma, mb)
mb
}
}
Main 用来输入两列数值与设置权重数据,并求得三次函数拟合的参数:
import scala.io.StdIn.readLine
object Main {
def main(args: Array[String]): Unit = {
var xc = Array[Double](0.00, 0.25, 0.50, 0.75, 1.00)
var yc = Array[Double](0.1, 0.35, 0.81, 1.09, 1.96)
var wc = Array[Double](1.0, 1.0, 1.0, 1.0, 1.0)
val cs = M2.sqrFit(xc, yc, wc)
cs.foreach(println(_))
readLine()
}
}
拟合结果如下,依次是常数、一次项参数、二次项参数、三次项参数: