Ok, so using Memoize.jl it’s not too bad until you get to the clincher… If you’re ok with a finite difference approximation, then it’s doable… let’s see what happens.
2 Likes
Using the hard stuff here, already;)
But, OTOH derivatives are somehow recursive too…
That’s what the finite differences approach leverages to get O(N) 
1 Like
Yes, it seems to be the way. I’ll will implement it and share the results here.
No need, start with this:
using Memoize
@memoize function wanddw(n,t,p)
if n < 0
error("wanddw: n must be an integer > 0")
end
if n == 0
if t < -1
return (complex(0.0),complex(0.0))
elseif t <= 1
return ((exp(p+p*t) - 1.0)/p, (t+1)*exp(p+p*t)/p - (exp(p+p*t)-1)/p^2)
else
return (2*exp(p*t)*sinh(p)/p, 2*sinh(p)*t*exp(p*t)/p - 2*sinh(p)*exp(p*t)/p^2 +2*cosh(p)*exp(p*t)/p)
end
elseif n == 1
t < -1 && return (complex(0.),complex(0.))
t <= 1 && return (-(1.0+exp(p+p*t)*(p-1)+p*t)/p^2, (-(p-1.0)*(t+1)*exp(p*t+p)-exp(p*t+p)-t)/p^2 - 2*((-(p-1)*exp(p*t+p)-p*t-1)/p^3))
return (-2.0*exp(p*t)*(p*cosh(p)-sinh(p))/p^2, -(2.0*(p*cosh(p)-sinh(p))*t*exp(p*t)/p^2)-2.0*sinh(p)*exp(p*t)/p - 4.0*(p*cosh(p)-sinh(p))*exp(p*t)/p^3)
elseif n >= 2
(wprev,dwprev) = wanddw(n-1,t,p)
(wprev2,dwprev2) = wanddw(n-2,t,p)
(wprev2d, dwprev2d) = wanddw(n-2,t,p+complex(0.0125,.0125))
return (2*t*wprev - 2*dwprev-wprev2, 2*t*dwprev - (dwprev2d-dwprev2)/complex(0.0125,0.0125))
end
end
The finite difference is hard coded, so you’ll want to fix that, but it returns IMMEDIATELY.
I used Maxima to get the symbolic derivatives for the base cases (because I’m lazy).
2 Likes
Thanks! I like the idea of returning both the value and the derivative in the same function for extra memoization efficiency!
I wonder if I can get this to work in terms of accuracy.
Thank you very much again.
It was fun! I think you just add the finite difference step size as a parameter and then use that in the code. Let us know how it works out, show some graphs or something.
I did that, and then tried it like this:
julia> wanddw(4,.1,complex(-.1),.0005)
(0.2959821757039954 - 0.0015914894716217987im, 0.4360000703131307 + 0.002182010051928633im)
julia> wanddw(4,.1,complex(-.1),.00005)
(0.2959821757039954 - 0.0015914894716217987im, 0.43592222276367454 + 0.0021033784685235563im)
julia> wanddw(4,.1,complex(-.1),.000005)
(0.2959821757039954 - 0.0015914894716217987im, 0.43591441315285184 + 0.002095691848207612im)
so at least as the step size gets smaller, it isn’t doing anything much
here’s the code with the dp:
using Memoize
@memoize function wanddw(n,t,p,dp)
if n < 0
error("wanddw: n must be an integer > 0")
end
if n == 0
if t < -1
return (complex(0.0),complex(0.0))
elseif t <= 1
return ((exp(p+p*t) - 1.0)/p, (t+1)*exp(p+p*t)/p - (exp(p+p*t)-1)/p^2)
else
return (2*exp(p*t)*sinh(p)/p, 2*sinh(p)*t*exp(p*t)/p - 2*sinh(p)*exp(p*t)/p^2 +2*cosh(p)*exp(p*t)/p)
end
elseif n == 1
t < -1 && return (complex(0.),complex(0.))
t <= 1 && return (-(1.0+exp(p+p*t)*(p-1)+p*t)/p^2, (-(p-1.0)*(t+1)*exp(p*t+p)-exp(p*t+p)-t)/p^2 - 2*((-(p-1)*exp(p*t+p)-p*t-1)/p^3))
return (-2.0*exp(p*t)*(p*cosh(p)-sinh(p))/p^2, -(2.0*(p*cosh(p)-sinh(p))*t*exp(p*t)/p^2)-2.0*sinh(p)*exp(p*t)/p - 4.0*(p*cosh(p)-sinh(p))*exp(p*t)/p^3)
elseif n >= 2
(wprev,dwprev) = wanddw(n-1,t,p)
(wprev2,dwprev2) = wanddw(n-2,t,p)
(wprev2d, dwprev2d) = wanddw(n-2,t,p+complex(dp,dp))
return (2*t*wprev - 2*dwprev-wprev2, 2*t*dwprev - (dwprev2d-dwprev2)/complex(dp,dp))
end
end
Even order 100 returns immediately, thanks to Memoization.jl !
julia> wanddw(100,.1,complex(-.1),.000005)
(6.329687336395586e86 - 2.2156952437460406e86im, -3.069266938528541e90 + 4.369343530297239e90im)
For what it’s worth, the above code is not numerically stable (and contains an error). I found another way to calculate the convolution analytically and use the recursive function to define the partial derivative with respect to p. Thanks again for the help!
1 Like