多项式插值用Julia怎么写,我这里出了点错误

我这里正准备超一遍 C 语言的多项式插值实现到Julia语言中,他的源代码在

这是我的 Julia 代码,他的运行有点错误

"""
  xs: 已知的点集\n
  fxs: xs 所对应的函数值\n
  zs: 所求的点集\n

  返回值: z 所对应的函数值
"""
function interpol(xs::Vector{T}, fxs::Vector{T}, zs::Vector{T}) where T <: Number
  coeff = Vector{Number}(undef, length(xs))
  table = Vector{Number}(undef, length(fxs))
  copyto!(table, fxs)
  pzs = Vector{Number}(undef, length(zs))
  n = length(xs)
  coeff[1] = table[1]
  
  for k in 2:n
    for i in 1:(n-k)
      j = i + k - 1
      table[i] = (table[i + 1] - table[i]) / (xs[j] - xs[i])
    end

    coeff[k] = table[1]
  end

  for k in 1:length(zs)
    pzs[k] = coeff[1]

    for j in 2:n
      term = coeff[j]
      for i in 1:j
        term = term * (zs[k] - xs[i])
      end

      pzs[k] = pzs[k] + term
    end
  end

  return pzs
end

他的测试代码是

@testset "test interpol" begin
  xs = 1:10
  fxs = map(x -> 1 + 2x + 3x^2, xs)
  @show fxs
  results = interpol(collect(xs), fxs, collect(11:20))
  @show results
  results = map(x -> 1 + 2x + 3x^2, 11:20)
  @show results
end

比较两个 results 的结果,他们应该相差不大才对,可是结果是


显然偏差极大
而C版本的就没有问题,测试代码是

int main() {
  double xs[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
  double fxs[] = {6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
  double zs[] = {11, 12, 13, 14, 15, 16, 17, 18, 19, 20};
  double pz[10];
  int result = interpol(xs, fxs, 10, zs, pz, 10);

  if (result == -1) {
    return -1;
  }

  for (int index = 0; index < 10; index += 1) {
    printf("%.2f\n", pz[index]);
  }
}

output is

gcc interpol.c && ./a.out 
386.00
457.00
534.00
617.00
706.00
801.00
902.00
1009.00
1122.00
1241.00

肯定是Julia的代码哪里错了,各位帮帮忙

哪里啊??。。。。。

你那边测试过了吗,测试过了再回复吧

做了点调试,似乎是第一次循环的问题
但是怎么看也看不出有什么问题,、

改为

      for i in 1:j-1
        term = term * (zs[k] - xs[i])
      end
1 个赞

是否调试过? 结果是否正确

你自己不会验证么?

卧槽,我验证完毕,正确的

作为求助者反而质问别人有没有验证过,得到正确结果连谢谢都不会说,还有这么不礼貌的人么?

5 个赞

备案号:京ICP备17009874号-2