一維二階切比雪夫多項式和二維二階切比雪夫多項式

scala版參考：

http://hxfcalf.blog.163.com/blog/static/21575548201373124214412

http://hxfcalf.blog.163.com/blog/static/21575548201373152715780

何老師的版本中的問題是採用了Float作為資料的儲存型別，但是由於Float的精度不夠，導致二維二階切比雪夫多項式結果不能收斂，

所以通過將資料型別修改為BigDecimal的方式解決了這個問題。

一維二階切比雪夫多項式：

package chebList
import tools.BigDecimalTools
/**
 
 * 一維2階切比雪夫多項式
 * gao
 * 2015-08-04
 */
case class CbVal(I0: Int) {
  val K0 = 2
//二階切比雪夫多項式係數
val cb = new Array[Seq[BigDecimal]](K0 + 1)
  for (k <- 0 to K0) cb(k) = k match{
    case 0 => (1 to I0).map (x => BigDecimal.decimal(1.0))
    case 1 => (1 to I0) map (x => BigDecimal.decimal(x - (I0 + 1.0
 
) / 2.0))
    case 2 => cb(1) map (x => x * x - (I0 * I0 - 1) / (4 * (4 - 1)))
  }

  //標準化處理後的二階切比雪夫多項式係數
val v = new Array[Seq[BigDecimal]](K0 + 1)
  def sq(v: Seq[BigDecimal]):BigDecimal = {
    BigDecimalTools.bigSqrt(v.map(x => x * x).sum)
  }
  for (k <- (0 to K0).par) {
    v(k) = cb(k).map(x => x / sq(cb
 
(k)))
  }
}

class D1Cheb(f: Array[BigDecimal]) extends CbVal(f.length) {
  //一維要素的多項式權重
val A = new Array[BigDecimal](K0 + 1)
  for (k <- (0 to K0).par) {
    A(k)= v(k).zip(f).map(x => x._1 * x._2).sum
  }
  //一維要素擬合
val fit = new Array[BigDecimal](f.length)
  for (i <- (0 until f.length).par){
    fit(i) = Seq(v(0)(i),v(1)(i),v(K0)(i)).zip(A).map(x => x._1 * x._2).sum
  }
  //一維要素擬合
def invers(len: Int, AA : Array[BigDecimal]):Array[BigDecimal]={
    val xx = new CbVal(len)
    val ffit = new Array[BigDecimal](len)
    for (i <- (0 until len).par){
      ffit(i) = Seq(xx.v(0)(i),xx.v(1)(i),xx.v(K0)(i)).zip(AA).map(x => x._1 * x._2).sum
    }
    ffit
  }
}

object CbVal extends App {
//  val c = CbVal(5)
//
//  for (k <- (0 to c.K0)) {
//    println(c.cb(k))
//  }
//  for (k <- (0 to c.K0)) {
//    println(c.v(k))
//  }
val data: Array[BigDecimal] = Array(0.0, 0.9999997, 2, 3.0, 4.0)
  val result:Array[BigDecimal] = Array(0.0,1.8,2.0,3.0,4.0)
  for(i <- (0 until 100)) {
    val curData = result
val cc = new D1Cheb(curData)
    val curResult = cc.fit
    result(1) = curResult(1)
    println(i,result(1))
  }
}

package chebList

/**
 * 二維2階切比雪夫多項式
 * gao
 * 2015-08-04
 */
import Array._
case class D2Cheb(f: Array[Array[BigDecimal]]) {
  val x = CbVal(f(0).length)
  val y = CbVal(f.length)
  //二維要素的多項式權重
val A = Array.fill(y.K0 + 1, x.K0 + 1)(BigDecimal.decimal(0.0))
  for (
    k <- (0 to y.K0); s <- (0 to x.K0);
i <- (0 until y.I0); j <- (0 until x.I0)
  ) {
    A(k)(s) += f(i)(j) * y.v(k)(i) * x.v(s)(j)
  }

  //二維要素擬合
val fit = Array.fill(y.I0, x.I0)(BigDecimal.decimal(0.0))
  for (
    i <- (0 until y.I0); j <- (0 until x.I0);
k <- (0 to y.K0); s <- (0 to x.K0)
  ) {
    fit(i)(j) += A(k)(s) * y.v(k)(i) * x.v(s)(j)
  }

  //二維要素擬合
def inverse(cols: Int, rows: Int, AA: Array[Array[BigDecimal]]): Array[Array[BigDecimal]] = {
    val xx = CbVal(cols)
    val yy = CbVal(rows)
    val z = ofDim[BigDecimal](yy.I0, xx.I0)
    for (
      i <- (0 until yy.I0).par; j <- (0 until xx.I0).par;
k <- (0 to yy.K0).par; s <- (0 to xx.K0).par
    ) {
      z(i)(j) += AA(k)(s) * yy.v(k)(i) * xx.v(s)(j)
    }
    z
  }
}

object D2Cheb extends App {
  val data: Array[Array[BigDecimal]] = Array(Array(0.0,1.0,2.0,3.0,4.0),
Array(1.0,2.6,3.0,4.0,5.0),
Array(2.0,3.0,4.0,5.0,6.0),
Array(3.0,4.0,5.0,6.0,7.0),
Array(4.0,5.0,6.0,7.0,8.0))

  for(i <- (0 until 100)){
    val curData = data
val cc = new D2Cheb(data)
    val curFit = cc.fit
    data(1)(1) = curFit(1)(1)
    println(i+", "+data(1)(1))
  }
}

python版本是基於scala版本進行改進的，能夠很好的收斂：

# -*- coding: utf-8 -*-
__author__ = 'Administrator'
'''
chebPolynomial
reference 周家斌，不規則格點上的車貝雪夫多項式展開問題 大氣科學 1983
'''
import numpy as np
import math

##二階切比雪夫多項式係數
def CbVal(IO):
    KO = 2
##二階切比雪夫多項式係數
cb = np.zeros([KO+1, IO])
    cb[0] = 1.0
cb1err = (IO+1.0)/2.0
cb[1] = [i-cb1err+1 for i in range(5)]
    cb2err = (IO*IO-1.0)/12.0
cb[2] = [i*i-cb2err for i in cb[1]]

    ##標準化處理後的二階切比雪夫多項式係數
v = np.zeros([KO+1, IO])
    sq0 = sq(cb[0])
    v[0] = [i/sq0 for i in cb[0]]
    sq1 = sq(cb[1])
    v[1] = [i/sq1 for i in cb[1]]
    sq2 = sq(cb[2])
    v[2] = [i/sq2 for i in cb[2]]
    return KO,cb,v

##求平方和的根
def sq(v):
    return math.sqrt(sum([i*i for i in v]))

##一維二階切比雪夫多項式
def D1Cheb(f):
    dimension = len(f)
    KO,cb,v = CbVal(dimension)

    A = np.zeros([KO+1])
    A[0] = sum(i[0]*i[1] for i in zip(f,v[0]))
    A[1] = sum(i[0]*i[1] for i in zip(f,v[1]))
    A[2] = sum(i[0]*i[1] for i in zip(f,v[2]))

    fit = np.zeros([dimension])
    for i in range(dimension):
        fit[i] = sum(i[0]*i[1] for i in zip([v[0][i],v[1][i],v[2][i]],A))
    return fit

##二維二階切比雪夫多項式
def D2Cheb(f):
    dimensionX = len(f[0])
    dimensionY = len(f)

    KOX,cbX,vX = CbVal(len(f[0]))
    KOY,cbY,vY = CbVal(len(f))

    ##二維要素的多項式權重
A = np.zeros([KOY+1, KOX+1])

    for k in range(KOY+1):
        for s in range(KOX+1):
            for i in range(dimensionY):
                for j in range(dimensionX):
                    A[k][s] += f[i][j] * vY[k][i] * vX[s][j]

    ##二維要素擬合
fit = np.zeros([dimensionY, dimensionX])
    for i in range(dimensionY):
        for j in range(dimensionX):
            for k in range(KOY+1):
                for s in range(KOX+1):
                    fit[i][j] += A[k][s] * vY[k][i] * vX[s][j]

    return fit

if __name__ == "__main__":
    # cb,v = CbVal(5)
    # print cb
    # print v
data = [[0.0, 1.0, 2.0, 3.0, 4.0],
[1.0, 2.6, 3.0, 4.0, 5.0],
[2.0, 3.0, 4.0, 5.0, 6.0],
[3.0, 4.0, 5.0, 6.0, 7.0],
[4.0, 5.0, 6.0, 7.0, 8.0]]

    fit = D2Cheb(data)

    for i in range(100):
        curData = data
        curFit = D2Cheb(curData)
        data[1][1] = curFit[1][1]
        print i, data[1][1]