我想搭建一个拉格朗日插值的算法程序来处理我的数据,锻炼我的算法能力。编程语言是 JuliaLang。
using DelimitedFiles
using Plots; pyplot()
function lagrange_interpolate(X,Y,t)
C = ones(length(X))
d = 0.0
for i = 1:length(X)
for j = [1:i-1;i+1:length(X)]
C[j] = C[j]*(t-X[j])/(X[i]-X[j])
end
d = d + Y[i] * C[i]
end
return d
end
A = readdlm("Numerical Methods/Data/data02.dat")
X = view(A,:,1)
Y = view(A,:,2)
T = 1.0:0.1:2.0
U = lagrange_interpolate.(X,Y,T)
plot([X;T],[Y;U])
savefig("U.png")
data02.dat:
0.0 0.0024979173609870897
0.1 0.03946950299855745
0.2 0.11757890635775581
0.3 0.22984884706593012
0.4 0.3662505856877064
0.5 0.5145997611506444
0.6 0.6616447834317517
0.7 0.7942505586276727
0.8 0.900571807773467
0.9 0.9711111703343291
1.0 0.9995675751366397
但是会得到错误的结果。 我想知道如何纠正它。
最佳答案
有两个问题。
首先,您的方法中存在错误。这是一个修复:
function lagrange_interpolate(X,Y,t)
C = ones(length(X))
d = 0.0
for i = 1:length(X)
for j = [1:i-1;i+1:length(X)]
C[i] = C[i]*(t-X[j])/(X[i]-X[j])
end
d = d + Y[i] * C[i]
end
return d
end
一个更简单的方法是写:
function lagrange_interpolate(X,Y,t)
idxs = eachindex(X)
sum(Y[i] * prod((t-X[j])/(X[i]-X[j]) for j in idxs if j != i) for i in idxs)
end
第二个问题是您错误地应用了广播。你应该写:
lagrange_interpolate.(Ref(X), Ref(Y), T)
因为您不希望X
和Y
被广播(并且Ref
保护包装在其中的值不被广播,参见https://docs.julialang.org/en/latest/manual/arrays/#Broadcasting-1 ).
在这种情况下也可能使用这样的理解:
[lagrange_interpolate(X, Y, t) for t in T]
会更容易阅读(但这是风格问题)。
关于algorithm - 如何更正我的 Julia 程序中的拉格朗日插值错误,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/55836257/