CODEWITHSUNDEEP
juliainterpolationautomatic-differentiation
HenryS
HenryS

Reputation: 21

Accurately finding the gradient of values from 2 vectors

So I have a vector of x-coordinates and a vector of y-coordinates and I need the second derivative at each x. There is no function tht I can write the y coordinates as therefore I don't think I can use forwarddiff.jl, finitdifferences.jl etc.

Thank you so much,

Henry

I have tried both BSplinekit.jl and Inteprolations.jl but this is linked to a long coupled ODE equation therefore the errors in the ends or centre of these methods seem to accumulate and after 1000 iterations lead to very eroneous results. Is there a possibilty that I could somehow use one of the function methods such as forwarddiff by having a package generate an 'exact' function from my results. If you need any data, please ask.

My result from Interpolations.jl after 2370-ish cycles (should be a straight line on 100) right now

enter image description here

My result from BSplineKit.jl after 1000-ish cycles

enter image description here

The code that I am using is as follows:

function Heatconductivity(κrt,Tel,Tph,zgrid)
    S=Interpolations.interpolate(zgrid,Tel,Gridded(Linear()))
    dS=only.(gradient.(Ref(S),zgrid))
    Q=dS.*kappa_linear(κrt,Tel,Tph)
    Qint=Interpolations.interpolate(zgrid,Q,Gridded(Linear()))
    dQ=only.(gradient.(Ref(Qint),zgrid))
    dQ=dQ.*-1
end

Where Tel is a vector of temperatures at depths into a material given by zgrid. I need the second derivative of the change in temperature with respect to position which currently I am getting from Interpolations.jl. I have tried multiple alogrithms for the interpolation but not seem to yield much improved results.

The output then is fed as part of a ODE that updates the temperature at each time step so I can't write an algebraic expression for dTel/dz.

Upvotes: 1

Views: 82

Answers (1)

Dan Getz
Dan Getz

Reputation: 18227

This is not the most accurate solution, but it may offer a quick stop-gap:

function d2f(x,y)
    h = diff(x)
    d2 = @. (2*(y[3:end] + h[2:end]/h[1:end-1]*y[1:end-2] - (1+h[2:end]/h[1:end-1])*y[2:end-1])/(h[2:end]*h[1:end-1]*(1+h[2:end]/h[1:end-1])))
    return d2
end

It is based on this old-ish paper: A simple finite-difference grid with non-constant intervals / Hilding Sundqvist & George Veronis (1970).

The returned vector is the second derivatives at points of x except for the endpoints (i.e. x[1] and x[end]).

It isn't a full interpolation method, but perhaps it is enough for your use-case.

Testing:

julia> x = [1.0:0.1:2.0...];

julia> y = x.^2;

julia> d2f(x,y)
9-element Vector{Float64}:
 2.000000000000005
 2.000000000000001
 2.000000000000005
 2.000000000000001
 2.0000000000000426
 1.9999999999999605
 2.0000000000000453
 1.9999999999999605
 1.9999999999999563

Indeed the second derivative of x^2 is 2.

Upvotes: 1

Related Questions