erf (Error function for complex numbers) - Visual Basic , VB.NET

	Find Code:
Home : Code : Forums : Submit : Mailing List : About : Contact

All

VB.NET

ASP.NET

C#

Snippets

Popular

Forums

Articles

Tips

Links

Books

Contest

erf (Error function for complex numbers)
Author: Paulo Buchsbaum
Website: http://www.greatsolutions.com.br/en
Submitted: 5/23/2012
Version: VB6
Compatibility: VB6
Category: Mathematics
Views: 12702
This function calculates erf (error function or gaussian aerror function) for complex and real numbers. erf ( x ) = 2 / sqrt (pi) * integration (t=0,x) [ exp ( - t ^2 ) dt ] There is no Visual Basic or VBA (Excel) code to calculate erf function for complex numbers, useful calculation for math applcations. I've based my code (see the source comments) in a code for MatLab, but it's more easy to do this in MatLab because MatLab does complex algebra. Gaussian Error Function is part of Cumulative Distribution Function (CDF) of a Normal Distribution: CDF = 1/2 * ( 1 + erf [ (x - m ) / raiz(2 s^2 ) ] Where m is the average, s is standard deviation and x is a aleatory variable. If a normal distribution has mean 0 and variance (s^2) = 0.5, the above expression becomes: CDF = 1/2 * (1 + erf(x) ) = 1/2 + 1/2 erf(x) The probability from -infinite to 0 is 1/2 from a normal distribution with mean 0. So erf(x) can be interpreted like twice the probability of a normal (average 0 and variance 0.5) to have a value between 0 and x.
Declarations:
'none
Code:
Option Explicit ************************************************************* '* Paulo Buchsbaum - May,23,2012 '* '* ERFZ Error function for complex (x,y). Return imaginary part. '* because for the most practical use in complex evaluation '* we use x=0 so the the real part is 0. '* The real part is returned on third parameter (RV) '* '* The function -erfz(0,i) solves the expression i * erf ( bi), because '* erf(bi) = 0 + ti, portanto i * erf(bi) =-ti '* '* .NET programmers - Includes Imports System.Math in you code '* Uses ATan and not Atn, and Sqrt and not Sqr '* (http://msdn.microsoft.com/en- us/library/thc0a116(v=vs.80).aspx) '* '* '* Based on MatLab Code (2001) written by '* Paul Godfrey <pgodfrey@intersil.com> and Peter J. Acklam <jacklam@math.uio.no> '* Ref: Abramowitz & Stegun section 7.1 equations 7.1.9, 7.1.23, and 7.1.29 '* (http://www.mathworks.com/matlabcentral/fileexchange/3574 -erfz) '************************************************** Function ErfZ(x As Double, y As Double, Optional ByRef Rv As Variant) As Double Const NMax = 193 Const nn = 32 Dim Re As Double Dim sqrtpi As Double Dim twopi As Double Dim az As Double Dim n As Integer Dim k1 As Double Dim k2r As Double, k2i As Double Dim k As Integer Dim xk As Double, yk As Double Dim Pi As Double Dim s1 As Double Dim s2r As Double, s2i As Double Dim s3 As Double Dim s4r As Double, s4i As Double Dim s5r As Double, s5i As Double Dim s6r As Double, s6i As Double Dim sr As Double, si As Double Dim yr As Double, yi As Double Dim a As Double, b As Double Dim Im As Double ' If real argument calculates with other function and exits If y = 0 Then ErfZ = ErfR(x) Exit Function End If Pi = 4 * Atn(1) ' Atn(1) = pi/4 (45 o.) twopi = 2 * Pi sqrtpi = Sqr(Pi) ' 1.772453850905516027298 az = Sqr(x ^ 2 + y ^ 2) If az <= 8 Then k1 = 2 / Pi * Exp(-x * x) ' Euler's Formula: e(xi) = cos(x) + i sin(x) ' ' k2 = Exp(-i * 2 * x * y) k2r = Cos(-2 * x * y) k2i = Sin(-2 * x * y) s1 = ErfR(x) If x <> 0 Then s2r = k1 / (4 * x) * (1 - k2r) s2i = k1 / (4 * x) * (1 - k2i) Else s2r = 0 s2i = 1 / Pi * y End If Re = s1 + s2r Im = s2i xk = x yk = y s5r = 0 s5i = 0 For n = 1 To nn s3 = Exp(-n * n / 4) / (n * n + 4 * xk * xk) 'Calculate s4 = 2xk - k2 (2xk cosh(nyk) - insinh(nyk)) ' ' (a+bi) * (c+di) = ' (ac-bd) + (ad + cb) i ' Where ' a = -k2r ' b = -k2i ' c = 2xk cosh(nyk) ' d = - n sinh(nyj) s4r = 2 xk - k2r * 2 * xk * Cosh(n * yk) - k2i * n * Sinh(n * yk) s4i = k2r * n * Sinh(n * yk) - k2i * 2 * xk * Cosh(n * yk) s5r = s5r + s3 * s4r s5i = s5i + s3 * s4i Next n s6r = k1 * s5r s6i = k1 * s5i Re = Re + s6r Im = Im + s6i Else If x <> 0 Then x = -x y = -y End If sr = 1 si = 0 ' Calculate y = 2 * z * z ' ' Complex squared ' ' ( a + bi ) * (a + bi ) = (a^2-b^2) + 2ab i yr = 2 * x ^ 2 - 2 * y ^ 2 yi = 2 * x * y For n = NMax To 1 Step -2 ' Complex Quotient ' ( a+ bi ) / (c + di) ' = (ac+bd)/(c2+d2) + (bc-ad)/(c2+d2)i ' ' Calculate s = 1 - n * (s / y) az = (yr ^ 2 + yi ^ 2) sr = 1 - n * (sr * yr + si * yi) / az si = -n * (si * yr - sr * yi) / az Next ' z = x + yi ' - z * z = - (x + yi) * (x+yi) = y^2 - x^2 - 2xy i ' ' Euler's Formula: e(xi) = cos(x) + i sin(x) ' ' Where exp(-z * z) = exp(y^2 - x^2) * exp(-2xyi) ' =exp(y^2 - x^2) * ( cos(-2xy) + i sin (-2xy) ) ' ' Complex Quotient ' ( a+ bi ) / (c + di) ' = (ac+bd)/(c2+d2) + (bc-ad)/(c2+d2)i ' ' Calculate f = 1 - s * Exp(-z * z) / (sqrtpi * z) ' ' Before calculate Exp(-z * z) / (sqrtpi * z) a = Exp(y ^ 2 - x ^ 2) * Cos(-2 * x * y) / sqrtpi b = Exp(y ^ 2 - x ^ 2) * Sin(-2 * x * y) / sqrtpi az = x ^ 2 + y ^ 2 Re = (a * x + b * y) / az Im = (b * x - a * y) / az ' Calculate now 1 -s * (Re,Im) ' ' Complex product: (a+bi) * (c+di) = (ac-bd) + (ad + bc) i ' Re = 1 - sr * Re + si * Im Im = -sr * Im - si * Re If x <> 0 Then Re = -Re Im = -Im End If If x = 0 Then Re = Re - 1 End If End If Rv = Re ErfZ = Im End Function '******************************************************* '* ErfReal '* Gaussian Error function with real argument X '* '* Based on '* Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables '* Milton Abramowitz e Irene A. Stegun '* Cited by John D. Cook that post a Python code in your blog '* ( http://www.johndcook.com/blog/2009/01/19/stand-alone- error-function-erf/ ) '********************************************************** Function ErfReal(x As Double) As Double ' Constants Const a1 = 0.254829592 Const a2 = -0.284496736 Const a3 = 1.421413741 Const a4 = -1.453152027 Const a5 = 1.061405429 Dim p As Double, t As Double, y As Double Dim sign As Integer p = 0.3275911 ' x signal sign = IIf(x < 0, -1, 1) x = Abs(x) t = 1 / (1 + p * x) y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Exp(-x * x) ErfReal = sign * y End Function '******************************************************* '* ErfR '* Gaussian Error function with real argument X '* using a standard Excel function '* '* Paulo Buchsbaum - May/2012 '******************************************************* Function ErfR(x As Double) As Double ' CDF of Normal Distribution = 1/2 * ( 1 + erf [ (x -m ) / raiz(2s^2 ) ] ' If Average = 0 and Standard Deviation = 1 ' CDF = 1/2 + 1/2 Erf (x) ' Erf = ( CDF-1/2 )2 (Double the probability from 0 to x ) ' On Error Resume Next ' Works on Excel's VBA ErfR = (WorksheetFunction.NormDist(x, 0, Sqr(1 / 2), True) - 0.5) * 2 If Err.Number <> 0 Then ErfR = ErfReal(x) On Error GoTo 0 End Function '******************************************************* '* Cosh '* Hyperbolic Cosine of x '***************************************************** Function Cosh(x As Double) As Double ' Modern VB Net has a native function math.cosh ' Excel users has worksheetfunction.cosh Cosh = (Exp(x) + Exp(-x)) / 2 End Function '***************************************************** '* Sinh '* Hyperbolic Syne of x '******************************************************* Function Sinh(x As Double) As Double ' Modern VB Net has a native function math.sinh ' Excel users has worksheetfunction.sinh Sinh = (Exp(x) - Exp(-x)) / 2 End Function