OpenOffice.org Forum at OOoForum.orgThe OpenOffice.org Forum
 
 [Home]   [FAQ]   [Search]   [Memberlist]   [Usergroups]   [Register
 [Profile]   [Log in to check your private messages]   [Log in

Kelvin Functions Ber(x), Bei(x), Ber'(x), Bei'(x)

 
Post new topic   Reply to topic    OOoForum.org Forum Index -> OpenOffice.org Code Snippets
View previous topic :: View next topic  
Author Message
UniqueUserID
General User
General User


Joined: 15 Sep 2009
Posts: 34

PostPosted: Mon Nov 02, 2009 12:52 am    Post subject: Kelvin Functions Ber(x), Bei(x), Ber'(x), Bei'(x) Reply with quote

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.
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    OOoForum.org Forum Index -> OpenOffice.org Code Snippets All times are GMT - 8 Hours
Page 1 of 1

 
Jump to:  
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


Powered by phpBB © 2001, 2005 phpBB Group