Как я могу исправить ошибку в моей программе Джулии для интерполяции Лагранжа - PullRequest
Новая
       56

Как я могу исправить ошибку в моей программе Джулии для интерполяции Лагранжа

1 голос
/ 24 апреля 2019

Я хочу построить программу алгоритма для интерполяции Лагранжа, чтобы обрабатывать мои данные и использовать мои способности алгоритма. Язык программирования - 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

Но это даст неверный результат. Я хочу знать, как это исправить.

My result

1 Ответ

2 голосов
/ 24 апреля 2019

Есть две проблемы.

Сначала у вас есть ошибка в вашем методе.Вот исправление: 

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]

будет легче читать (но это вопрос стиля).

...