我这里正准备超一遍 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的代码哪里错了，各位帮帮忙