Defining a function using the output of another function

General Usage
, ,
1

I have a function that takes an array and returns its derivative.

function dy(y, δx)

    N = length(y)
    dy = zeros(N)
    for j = 4:N-3
        dy[j]=( 5*(y[j+1] - y[j-1])
            + 4*(y[j+2] - y[j-2])
            + 1*(y[j+3] - y[j-3]))/(32*δx)
    end

    j=3
    dy[j]=( 2*(y[j+1] - y[j-1])
        + 1*(y[j+2] - y[j-2]))/(8*δx)

    j=N-2
    dy[j]=( 2*(y[j+1] - y[j-1])
        + 1*(y[j+2] - y[j-2]))/(8*δx)

    j=2
    dy[j]=(y[j+1] - y[j-1])/(2*δx)

    j=N-1
    dy[j]=(y[j+1] - y[j-1])/(2*δx)

    dy[1]=2*dy[2]-dy[3]
    dy[N]=(y[N-2]-4*y[N-1]+3*y[N])/(2*δx)

    return dy

end

There are basically two things required to define such a function:

  1. The coefficients multiplying the values [-1,-4,-5,5,4,1]./12
  2. The corresponding array positions wrt to the current element [-3,-2,-1,1,2,3]

I have written a different function which gives me these two things for a given order of derivative and stencil size. For example, it will return the following arrays for second-order derivate using five stencil points

Coefficients
2.91667 -8.66667 9.5 -4.66667 0.916667
0.916667 -1.66667 0.5 0.333333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 1.33333 -2.5 1.33333 -0.0833333
-0.0833333 0.333333 0.5 -1.66667 0.916667
0.916667 -4.66667 9.5 -8.66667 2.91667

Relative positions
0 1 2 3 4
-1 0 1 2 3
-2 -1 0 1 2
-2 -1 0 1 2
-2 -1 0 1 2
-2 -1 0 1 2
-2 -1 0 1 2
-2 -1 0 1 2
-2 -1 0 1 2
-3 -2 -1 0 1
-4 -3 -2 -1 0

Is there a way to use this information to write the function dy ?

One dirty way is to simply use the above arrays and string manipulation to write a file with the new version of dy during execution and then include it.
But, I think there should be a better way to do this.

2

I’m not sure if I understood the problem correctly, but perhaps just pass the coefficients and offsets as parameters? Something like

function dy(y, δx, coeffs, offsets)
    N = length(y)
    dy = zeros(N)
    for j in eachindex(y)
        denom = 32*δx  # and something else near the edges?
        dy[j] = sum(y[j+i]*c for (i,c) in zip(offsets,coeffs) if j+i in eachindex(y)) / denom
    end
    dy
end

(I’m not sure what that denominator is for and where the numbers 32,8,2 come from).
You might want to take a look at ImageFiltering.jl, and specifically ImageFiltering.imfilter. It does exactly what it seems you are trying to do, except maybe near the edges (since I’m not sure what is going on near the edges in your code).