Skip to content

Commit 6c12c23

Browse files
committed
add taylor series
1 parent edbd576 commit 6c12c23

File tree

1 file changed

+50
-0
lines changed

1 file changed

+50
-0
lines changed

src/besselj.jl

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -187,3 +187,53 @@ function besselj(n::Int, x::Float64)
187187

188188
return sign * ans
189189
end
190+
191+
192+
function b2(v, x)
193+
194+
a = (x/2)^v
195+
out = 1 / gamma(v + 1)
196+
for k in 1:100
197+
out += (-x^2 / 4)^k / factorial(BigFloat(k)) / factorial(v + BigFloat(k))
198+
end
199+
return out*a
200+
end
201+
202+
function b(nu, x::T) where T
203+
coef = (x/2)^nu / gamma(nu + 1)
204+
205+
x2 = (x/2)^2
206+
207+
208+
maxiter = 10000
209+
b = one(T)
210+
for k in 1:maxiter
211+
z = -x2 / (k * (k + nu))
212+
b += z*b
213+
end
214+
return b * coef
215+
end
216+
217+
218+
### for this function dd can be calculated accurately but there is a build up in error for val as we sum
219+
## looks like only works for small arguments and smallish orders....
220+
function b3(v, x::T) where T
221+
MaxIter = 5000
222+
if v < 20
223+
coef = (x / 2)^v / factorial(v)
224+
else
225+
vinv = inv(v)
226+
coef = sqrt(vinv / (2 * π)) * MathConstants.e^(v * (log(x / (2 * v)) + 1))
227+
coef *= evalpoly(vinv, (1, -1/12, 1/288, 139/51840, -571/2488320, -163879/209018880, 5246819/75246796800, 534703531/902961561600))
228+
end
229+
230+
x2 = (x / 2)^2
231+
val = zero(T)
232+
233+
for k in 1:MaxIter
234+
val += coef
235+
coef = -coef * x2 / (k * (v + k))
236+
abs(coef) < eps(T) * abs(val) && break
237+
end
238+
return val
239+
end

0 commit comments

Comments
 (0)