[Home]   [FAQ]   [Search]   [Memberlist]   [Usergroups]   [Register]

Author Message
UniqueUserID
General User

Joined: 15 Sep 2009
Posts: 34

Posted: Mon Nov 02, 2009 12:52 am    Post subject: Kelvin Functions Ber(x), Bei(x), Ber'(x), Bei'(x)

If you ever need these functions, this should save you a bit of time:

 Code: REM ****** Kelvin Functions ********* Function Ber(ByVal x as double) as double   ' Calculates Ber(x) function   ' Uses recurrence relation based on series expansion formula 820.3   ' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"   ' 4th Edition, 1961, MacMillan    ' Basic code Copyright 2009, Robert Weaver   ' Revised: 2009-11-13   ' Permission is hereby granted to copy and distribute verbatim copies of this software.   ' This software is without warranty. The user accepts all risks incurred from its use.   Dim i As integer,sign As integer   Dim sum As Double,termi As Double,term As Double   if x=0 then     Ber()= 1   else     term=1     sign=1     sum=term*sign     '300 iterations is enough to calculate any value within the range of double precision     for i = 1 to 300       sign=-sign       termi=((x/(4*i)) * (x/(4*i-2)))^2       term=term*termi       sum=sum+sign*term       'Skip out of loop if current term is less than 1e-12 of value       if abs(term/sum)<1e-12 then exit for     next     Ber()= sum   end if end Function Function Ber_(ByVal x as double) as double   ' Calculates Ber'(x) function -- d/dx Ber(x)   ' Uses recurrence relation based on series expansion formula 820.5   ' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"   ' 4th Edition, 1961, MacMillan   ' Basic code Copyright 2009, Robert Weaver   ' Revised: 2009-11-13   ' Permission is hereby granted to copy and distribute verbatim copies of this software.   ' This software is without warranty. The user accepts all risks incurred from its use.   Dim i  As integer,sign As integer   Dim sum As Double,termi As Double,term As Double   if x=0 then     Ber_()= 0   else     term=(x/2)^3/2     sign=-1     sum=term*sign     '300 iterations is enough to calculate any value within the range of double precision     for i = 1 to 300       sign=-sign       termi=(x/(4*i)) * (x/(4*i+2))^2 * (x/(4*i+4))       term=term*termi       sum=sum+sign*term       'Skip out of loop if current term is less than 1e-12 of value       if abs(term/sum)<1e-12 then exit for     next     Ber_()= sum   end if end Function Function Bei(ByVal x as double) as double   Dim i  As integer,sign As integer   Dim sum As Double,termi As Double,term As Double   ' Calculates Bei(x) function   ' Uses recurrence relation based on series expansion formula 820.4   ' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"   ' 4th Edition, 1961, MacMillan   ' Basic code Copyright 2009, Robert Weaver   ' Revised: 2009-11-13   ' Permission is hereby granted to copy and distribute verbatim copies of this software.   ' This software is without warranty. The user accepts all risks incurred from its use.   if x=0 then     Bei()= 0   else     term=x*x/4     sign=1     sum=term*sign     '300 iterations is enough to calculate any value within the range of double precision     for i = 1 to 300       sign=-sign       termi=((x/(4*i)) * (x/(4*i+2)))^2       term=term*termi       sum=sum+sign*term       'Skip out of loop if current term is less than 1e-12 of value       if abs(term/sum)<1e-12 then exit for     next      Bei()= sum   end if end Function Function Bei_(ByVal x as double) as double   Dim i  As integer,sign As integer   Dim sum As Double,term As Doublei,term As Double   ' Calculates Bei'(x) function -- d/dx Bei(x)   ' Uses recurrence relation based on series expansion formula 820.6   ' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"   ' 4th Edition, 1961, MacMillan   ' Basic code Copyright 2009, Robert Weaver   ' Revised: 2009-11-13   ' Permission is hereby granted to copy and distribute verbatim copies of this software.   ' This software is without warranty. The user accepts all risks incurred from its use.   if x=0 then      Bei_()= 0   else     term=x/2     sign=1     sum=term*sign     '300 iterations is enough to calculate any value within the range of double precision     for i = 1 to 300       sign=-sign       termi=(x/(4*i-2)) * (x/(4*i))^2 * (x/(4*i+2))       term=term*termi       sum=sum+sign*term        'Skip out of loop if current term is less than 1e-12 of value       if abs(term/sum)<1e-12 then exit for     next      Bei_()= sum   end if end Function

2009-11-13, Edit Note:
I recently discovered that my original dimension statements

Dim i,sign As integer
Dim sum,termi,term As Double

do not declare ALL variables on the respective line to be all integers or all doubles, just the variable immediately preceding the 'As'.

Therefore I've edited the Dim statements so that all variables are declared correctly.

The failure to declare the floating point variables as doubles didn't seem to have hurt accuracy however. I compared the results of original versions of the functions to published tables and they agree to 10 decimal places. However, it's safer to have things declared properly.
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT - 8 Hours Page 1 of 1

 Jump to: Select a forum OpenOffice.org Forums----------------Setup and TroubleshootingOpenOffice.org WriterOpenOffice.org CalcOpenOffice.org ImpressOpenOffice.org DrawOpenOffice.org MathOpenOffice.org BaseOpenOffice.org Macros and APIOpenOffice.org Code Snippets Community Forums----------------General DiscussionSite Feedback
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum