二阶切比雪夫多项式实现（scala版、python版）

一维二阶切比雪夫多项式和二维二阶切比雪夫多项式

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]