View Full Version : how do you solve a system of linear equations in ThinBasic?
in the thread by Dan on the Gamma function, there's mention of the Lanczos approximation.
the Lanczos coefficients can be calculated by solving a system of linear equations but I need
help on how to do it in ThinBasic.
in the example given in Numerical Recipies, n = 7, g = 5
so you set up a matrix like this
'for n = 7, the following sets up matrix a as:
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=1 To n
a(i,1) = 1
For j=2 To n
a(i,j)=1/(i+j-2)
Next
Next
then you set up the equation vector like this
For i=1 To n
c(i) = 0.5*Factorial(i-1)/((i+g-0.5)^(i-0.5))*Exp(i+g-0.5)*Sqr(2)/Sqr(M_PI)
Next
and you find the coefficients by solving a() = c()
in Maple I would do it like this
multiply(inverse(a),c)
here's my attempt in ThinBasic, for more terms and accuracy you need much higher precision.
Uses "console"
Uses "math"
Const n = 7
Dim As Ext g = 5
Dim As Ext a(n,n), b(n,n), c(n,1)
Dim As Integer i, j
'for n = 7, the following sets up matrix a as:
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=1 To n
a(i,1) = 1
For j=2 To n
a(i,j)=1/(i+j-2)
Next
Next
For i=1 To n
c(i,1) = Factorial(i-1)/((i+g-0.5)^(i-0.5)/Exp(i+g-0.5)*Sqr(2*M_PI))
Next
MAT b() = INV(a())
MAT a() = b() * c()
For i=1 To n
Console_WriteLine "C("&i&") = "&a(i,1)
Next
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
admittedly you will need arbitrary precision to get good results.
danbaron
13-08-2011, 07:08
I don't know where you got the values for matrix a from.
Anyway, I'm getting that it is singular, in which case, as you know, it has no inverse.
Dan
' code --------------------------------------------------------------------------------
Uses "console"
Uses "math"
Const n = 7
Dim As Ext a(n,n), b(n,n)
Dim As Integer i, j
Dim As Ext dta, dtb
'for n = 7, the following sets up matrix a as:
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=1 To n
a(i,1) = 1
For j=2 To n
a(i,j)=1/(i+j-2)
Next
Next
' Check the DET function.
randomize
For i = 1 To n
For j = 1 To n
b(i,j) = Rnd * (i^2 + j^2)
Next
Next
dta = DET(a())
dtb = DET(b())
Console_WriteLine "The determinant of a = " & dta
Console_WriteLine "The determinant of b = " & dtb
Console_WriteLine
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
' output ------------------------------------------------------------------------------
The determinant of a = 5.80875406944508548E-21
The determinant of b = 156682697.391266545
All done. Press any key to finish
Dan, I beg to disagree as to the matrix a being singular, ThinBasic will calculate the inverse OK though not very precise,
and Maple solves it without a problem, but I am unsure whether the row and column order are the same in ThinBasic as in Maple.
as to where I got the numbers for the matrix frankly I forgot, it's been like 20 year since I figured out how to set up the
equations.
here's how you calculate the Lanczos coefficients in Maple, first compute the coefficients as given in Numerical Recipies
note: I put a Maple comment marker on the Maple output.
> restart;
> Digits := 40;
> with(linalg);
> f1 := proc (n) local i, j, v; v := matrix(n+1, n+1); for i to n+1 do v[i, 1] := 1; for j from 2 to n+1 do v[i, j] := 1/(i+j-2) end do end do; RETURN(v) end proc;
> z := unapply(factorial(i-1)*exp(i+g-1/2)/((i+g-1/2)^(i-1/2)*sqrt(2*Pi)), i, g);
> f2 := proc (n, y) local i, v; v := matrix(n+1, 1); for i to n+1 do v[i, 1] := z(i, y) end do; RETURN(v) end proc;
> n := 7;
> g := 5;
> f := f1(n);
> g1 := inverse(f);
> b := f2(n, g);
> gg := multiply(g1, b);
> evalf(gg[1, 1]-1, 4*rhs(op(2, op(1, gg))[1])+4);
#output -9.3304068958898243168 e-12
> evalf(gg[rhs(op(2, op(1, gg))[1]), 1], 4*rhs(op(2, op(1, gg))[1])+4);
#output 0.000002394534913406850300243
> for i to rhs(op(2, op(1, gg))[1]) do evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+2*i) end do;
#output
# 0.99999999999066959310411017566
# 76.1800917308668800992939735513657
# -86.5053203963967652303139216042357
# 24.01409899435588324058704081692617
# -1.231742921450034278218762025469663
# 0.00121555828611639057422612294568442
# -0.0000120262591451567138253281149246451
# 0.000002394534913406850300240770820124155
> f3 := evalf(gg[1, 1], 4*rhs(op(2, op(1, gg))[1])+1);
> for i from 2 to rhs(op(2, op(1, gg))[1]) do f3 := f3+evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+i)/(x+i-1) end do;
> f4 := unapply((x+g+1/2)^(x+1/2)*exp(-x-g-1/2)*sqrt(2*Pi)*f3, x);
> evalf(f4(1/2));
#output
# 0.8862269254527589674644560128649165283764
> 4*(%*%)-evalf(Pi);
#output
# 6.762374918539830701792077 e-15
now compute the coefficients I cited.
> restart;
> Digits := 40;
> with(linalg);
> f1 := proc (n) local i, j, v; v := matrix(n+1, n+1); for i to n+1 do v[i, 1] := 1; for j from 2 to n+1 do v[i, j] := 1/(i+j-2) end do end do; RETURN(v) end proc;
> z := unapply(factorial(i-1)*exp(i+g-1/2)/((i+g-1/2)^(i-1/2)*sqrt(2*Pi)), i, g);
> f2 := proc (n, y) local i, v; v := matrix(n+1, 1); for i to n+1 do v[i, 1] := z(i, y) end do; RETURN(v) end proc;
> n := 14;
> g := 607/128;
> f := f1(n);
> g1 := inverse(f);
> b := f2(n, g);
> gg := multiply(g1, b);
> evalf(gg[1, 1]-1, 4*rhs(op(2, op(1, gg))[1])+4);
#output
# -2.90817953577301957618515490814322962380380 e-15
> evalf(gg[rhs(op(2, op(1, gg))[1]), 1], 4*rhs(op(2, op(1, gg))[1])+4);
#output
# 0.00000368991826595316227036759674456787362014192
> for i to rhs(op(2, op(1, gg))[1]) do evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+2*i) end do;
#output
# 0.999999999999997091820464226980423814845091856770376328
# 57.156235665862923516579393485977444273254442892929107852
# -59.597960355475491248142266131602836001910100856036329460
# 14.1360979747417471738634195408767260156591641388059836898
# -0.4919138160976201997828400285303078441634827850503476487
# 0.00003399464998481188869891934155234155285736767075962214
# 0.000046523628927048575665230224965822889358967609346662905
# -0.00009837447530487956467653837063470730301913302874861670310
# 0.000158088703224912488836072413443663958858739583487234932827
# -0.00021026444172410488319269928300571190544166498661591073776001
# 0.0002174396181152126431961446496038346923763748760822190719415235
# -0.000164318106536763890217069562280184132420847097097457818988636042
# 0.00008441822398385274329281181534545498133890999303130609738370698039
# -0.0000261908384015814086696650362479549018414325910258739195513825336890
# 0.0000036899182659531622703675967445678736201430201881942403488492031328990
> f3 := evalf(gg[1, 1], 4*rhs(op(2, op(1, gg))[1])+1);
> for i from 2 to rhs(op(2, op(1, gg))[1]) do f3 := f3+evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+i)/(x+i-1) end do;
> f4 := unapply((x+g+1/2)^(x+1/2)*exp(-x-g-1/2)*sqrt(2*Pi)*f3, x);
> evalf(f4(1/2));
#output
# 0.8862269254527580134320931844229516344461
# the % character is shorthand for last result
> 4*(%*%)-evalf(Pi);
#output
# -1.538422995214718383218 e-18
danbaron
13-08-2011, 20:36
I don't have much time now, but, I confirmed you are right, Johan.
The inverse of a is correct.
Dan
' code --------------------------------------------------------------------------------
Uses "console"
Uses "math"
Const n = 7
Dim As Ext a(n,n), b(n,n), c(n,n)
Dim As Integer i, j
Dim As Ext dta, dtb
'for n = 7, the following sets up matrix a as:
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=1 To n
a(i,1) = 1
For j=2 To n
a(i,j)=1/(i+j-2)
Next
Next
MAT b() = INV(a())
MAT c() = a() * b()
For i = 1 To n
For j = 1 To n
PrintL c(i,j)
Next
PrintL
Next
Console_WriteLine
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
' output ------------------------------------------------------------------------------
1.00000000000000039
-7.10542735760100186E-15
5.68434188608080149E-14
-4.54747350886464119E-13
4.54747350886464119E-13
-4.54747350886464119E-13
5.68434188608080149E-14
2.22044604925031308E-16
.999999999999994671
5.68434188608080149E-14
-3.97903932025656104E-13
3.41060513164848089E-13
-1.1368683772161603E-13
8.52651282912120223E-14
2.77555756156289135E-16
-5.32907051820075139E-15
1.00000000000005684
-2.84217094304040074E-13
3.41060513164848089E-13
-5.68434188608080149E-14
1.1368683772161603E-13
1.11022302462515654E-16
0
1.42108547152020037E-14
.999999999999772626
2.27373675443232059E-13
-1.1368683772161603E-13
5.68434188608080149E-14
2.22044604925031308E-16
-3.55271367880050093E-15
5.68434188608080149E-14
-2.27373675443232059E-13
1.00000000000022737
-5.68434188608080149E-14
2.84217094304040074E-14
1.38777878078144568E-16
-3.55271367880050093E-15
1.42108547152020037E-14
-1.1368683772161603E-13
0
.999999999999943157
2.84217094304040074E-14
8.32667268468867405E-17
0
-1.42108547152020037E-14
-5.68434188608080149E-14
1.1368683772161603E-13
1.1368683772161603E-13
1.00000000000002842
All done. Press any key to finish
Petr Schreiber
13-08-2011, 21:05
Hi guys,
I have one testing example of linear equations solution on my drive. The following problem:
1x1 + 1x2 + 1x3 = 23
0x1 + 1x2 - 1x3 = -1
1x1 + 0x2 - 2x3 = 0
can be solved using this code, please note the use of square brackets to make possible write matrices in row order:
Uses "Math"
Dim A(3, 3) As Ext
Dim b(3, 1) As Ext
Dim x(3, 1) As Ext
A(1,1) = [ 1, 1, 1,
0, 1,-1,
1, 0,-2 ]
b(1,1) = [ 23,
-1,
0 ]
MAT A() = INV( A() )
MAT x() = A() * b()
MsgBox 0, x(1, 1)
MsgBox 0, x(2, 1)
MsgBox 0, x(3, 1)
Petr
thanks Petr
in my example the matrix is such that you need high precision, but it's good to know how to solve in ThinBasic.
:)
fixed my first post.
REDEBOLT
14-08-2011, 06:49
Hi Jack,
Allow me to have a go at this.
Here is the solution as calculated by Octave-3.2.4 ( a MATLAB clone ).
octave-3.2.4.exe:2> A=[
> 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6
> 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7
> 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8
> 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9
> 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10
> 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11
> 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12
> ]
A =
1.000000 1.000000 0.500000 0.333333 0.250000 0.200000 0.166667
1.000000 0.500000 0.333333 0.250000 0.200000 0.166667 0.142857
1.000000 0.333333 0.250000 0.200000 0.166667 0.142857 0.125000
1.000000 0.250000 0.200000 0.166667 0.142857 0.125000 0.111111
1.000000 0.200000 0.166667 0.142857 0.125000 0.111111 0.100000
1.000000 0.166667 0.142857 0.125000 0.111111 0.100000 0.090909
1.000000 0.142857 0.125000 0.111111 0.100000 0.090909 0.083333
octave-3.2.4.exe:3> format long g
octave-3.2.4.exe:4> det(A)
ans = 5.80876610909398e-021
octave-3.2.4.exe:5> B=inv(A)
B =
Columns 1 through 3:
1.00000000090092 -42.0000000317041 420.000000271976
42.0000000068425 -882.000000243734 5880.0000021059
-840.000000232474 23520.0000082602 -176400.000071267
5040.00000181503 -158760.000064399 1270080.0005551
-12600.0000053416 423360.000189264 -3528000.00163032
13860.0000065701 -485100.000232675 4158000.00200328
-5544.00000285337 199584.000100984 -1746360.00086911
Columns 4 through 6:
-1680.00000095145 3150.00000158441 -2772.00000125375
-17640.0000074021 26460.0000123685 -19404.0000098124
564480.000250259 -882000.000417884 665280.000331352
-4233600.00194806 6804000.00325144 -5239080.00257731
12096000.005719 -19845000.0095424 15523200.0075621
-14553000.007025 24255000.0117187 -19209960.0092852
6209280.00304696 -10478160.0050819 8382528.00402599
Column 7:
924.000000379607
5544.00000297698
-194040.000100488
1552320.00078141
-4656960.00229232
5821200.00281424
-2561328.0012201
octave-3.2.4.exe:6> C=[
> 41.6244369164390698
> 16.012316405251682
> 9.36473553710404886
> 6.57049625931606175
> 5.095216729296468
> 4.20380175313001118
> 3.61487381446333504
> ]
C =
41.6244369164391
16.0123164052517
9.36473553710405
6.57049625931606
5.09521672929647
4.20380175313001
3.61487381446334
octave-3.2.4.exe:7> D=B*C
D =
1.0000000002019
76.1800917295827
-86.5053203336004
24.0140982591429
-1.23173967863841
0.00120870569662657
-5.41476765647531e-006
octave-3.2.4.exe:8> diary off
And here is the program:
' Empty GUI script created on 08-13-2011 03:38:47 by Robert E. DeBolt (ThinAIR)
Uses "console"
Uses "math"
Const n = 7
Dim As Ext g = 5
Dim As Ext a(n,n), b(n,n), c(n,1), d(n,1)
Dim As Integer i, j
'for n = 7, the following sets up matrix a as:
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=1 To n
a(i,1) = 1
For j=2 To n
a(i,j)=1/(i+j-2)
Next
Next
For i=1 To n
c(i,1) = 0.5*Factorial(i-1)/((i+g-0.5)^(i-0.5))*Exp(i+g-0.5)*Sqr(2)/Sqr(M_PI)
PrintL c(i,1)
Next
PrintL
MAT b() = INV(a())
For i=1 To n
For j=1 To n
' Print "b(" & i & Str$(j) & ") = " & str$(b(i,j),5) & " "
Print Str$(b(i,j),10) & $TAB
Next j
PrintL
Next i
MAT d() = b() * c()
For i=1 To n
Console_WriteLine "D("& i &") = "&d(i,1)
Next i
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
And here are the results:
41.6244369164390698
16.012316405251682
9.36473553710404886
6.57049625931606175
5.095216729296468
4.20380175313001118
3.61487381446333504
1 -42 420 -1680 3150 -2772 924
42 -882 5880 -17640 26460 -19404 5544
-840 23520 -176400 564480 -882000 665280 -194040
5040 -158760 1270080 -4233600 6804000 -5239080
1552320
-12600 423360 -3528000 12096000 -19845000 15523200
-4656960
13860 -485100 4158000 -14553000 24255000 -19209960
5821200
-5544 199584 -1746360 6209280 -10478160 8382528
-2561328
D(1) = 1.00000000019002488
D(2) = 76.1800917294715312
D(3) = -86.505320329419078
D(4) = 24.0140982408456694
D(5) = -1.23173957250764943
D(6) = 1.20865104327094741E-3
D(7) = -5.39526990905869752E-6
All done. Press any key to finish
Regards,
Bob
danbaron
14-08-2011, 07:57
(I didn't see what Bob posted until I finished.)
Anyway, maybe you changed something.
Because, I thought that before when I ran it, it gave crazy answers.
That's why I checked to see if the determinant of matrix "a" was 0.
Now, however, when I run it, it seems OK (maybe I did something wrong before). Below, I'll compare the output, with the coefficients from "Numerical Recipes in C", 2nd ed., p. 214, function gammaln.
(Numerical Recipes does not calculate the coefficient with value 1.)
' code ----------------------------------------------------------------------------------
Uses "console"
Uses "math"
Const n = 7
Dim As Ext g = 5
Dim As Ext a(n,n), b(n,n), c(n,1)
Dim As Integer i, j
'for n = 7, the following sets up matrix a as:
' [[ 1, 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]
' [ 1, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12 ]]
For i=1 To n
a(i,1) = 1
For j=2 To n
a(i,j)=1/(i+j-2)
Next
Next
For i=1 To n
c(i,1) = Factorial(i-1)/((i+g-0.5)^(i-0.5)/Exp(i+g-0.5)*Sqr(2*M_PI))
Next
MAT b() = INV(a())
MAT a() = b() * c()
For i=1 To n
Console_WriteLine "C("&i&") = "&a(i,1)
Next
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
' output --------------------------------------------------------------------------------
C(1) = 1.00000000019002555
C(2) = 76.1800917294715294
C(3) = -86.5053203294193054
C(4) = 24.0140982408470336
C(5) = -1.2317395725131064
C(6) = 1.20865105054690503E-3
C(7) = -5.39527172804810107E-6
All done. Press any key to finish
' Numerical Recipes coefficients --------------------------------------------------------
cof[0] = 76.18009172947146;
cof[1] = -86.50532032941677;
cof[2] = 24.01409824083091;
cof[3] = -1.231739572450155;
cof[4] = 0.1208650973866179e-2;
cof[5] = -0.5395239384953e-5;
Check the differences between the two versions of the coefficients.
' code ----------------------------------------------------------------------------------
Uses "console"
Uses "math"
Const n = 6
Dim As Ext a(n), b(n), c(n)
Dim As Integer i
' From Johan.
a(1) = 76.1800917294715294
a(2) = -86.5053203294193054
a(3) = 24.0140982408470336
a(4) = -1.2317395725131064
a(5) = 1.20865105054690503 * 10^-3
a(6) = -5.39527172804810107 * 10^-6
' From Numerical Recipes.
b(1) = 76.18009172947146
b(2) = -86.50532032941677
b(3) = 24.01409824083091
b(4) = -1.231739572450155
b(5) = 0.1208650973866179 * 10^-2
b(6) = -0.5395239384953 * 10^-5
' abs(difference).
For i = 1 To n
c(i) = Abs(a(i)-b(i))
Next
For i = 1 To n
PrintL c(i)
Next
Console_WaitKey
' output --------------------------------------------------------------------------------
6.9395877932976191E-14
2.53540244354866218E-12
1.61235989837305027E-11
6.29514000606221091E-11
7.66807260298555909E-11
3.23430951010704571E-11
In thinBasic, pi is defined to 16 digits.
I tried instead defining it myself,
dim as ext pie = 4*atn(1)
but, it didn't help.
Type EXT is supposed to have 18 digits of precision, but, apparently your calculations are not going to come very close to that.
The one thing I can think of to try in order to get more accuracy, is to solve the system by Gaussian Elimination instead of by inverting "a". We can't see how much accuracy is lost by using the inversion function, but we can get an idea. If you look in my previous post, you see that a * INV(a) in the worst cases is only accurate to 12 digits.
On the other hand, now that I think about it, the lack of accuracy in the calculation of INV(a), is probably just be because "a" is ill-conditioned. When I calculated its determinant, I got,
The determinant of a = 5.80875406944508548E-21
I think that number is pretty bad for the determinant of a matrix. A matrix which cannot be inverted has determinant 0. The closer the determinant of the coefficient matrix is to 0, the worse will be the results when the coefficient matrix (or its inverse) is used to solve a system of linear equations.
REDEBOLT
14-08-2011, 11:14
Hi Dan,
The determinant calculation from Octave gives
ans = 5.80876610909398e-021 from Octave
which agrees with the Thinbasic calculation to about five significant digits.
The coefficient matrix is called a Hilbert (http://en.wikipedia.org/wiki/Hilbert_matrix) matrix which is notoriously ill-conditioned, but is of great theoretical interest in the study of gamma functions.
I have big integer package which can calculate with arbitrary precision. I will see what I can do with the Gauss-Jordan elimination procedure for calculating the matrix inverse.
Regards,
Bob
danbaron
14-08-2011, 11:35
Well, good then, Bob.
I know of David Hilbert.
I'm happy to have your confirmation concerning the ill-conditioning.
Dan
----------------------------------------------------------------
Now, I looked at your reference.
I like it -->, "making them notoriously difficult to use in numerical computation".
----------------------------------------------------------------
I think that no matter how ill-conditioned a matrix is, i.e., no matter how close the determinant is to 0, it is still theoretically possible to get good numerical results when using it to solve a system of linear equations.
But, in order to do so, you may need a million digits of precision.
I think the only case that is hopeless, is when the matrix is singular, the determinant is exactly 0. Then, the n x n matrix does not contain enough information to uniquely solve for n unknowns.
hi Bob and Dan
Bob, Octave does much beter than thinBasic, Dan, I did change the script to make it work.
anyway, you can also drop the first column in the matrix and just make the first coefficient = 1,
(also decrease the size of the dimensions by one) the rest of the coefficients will differ
slightly because of this change, but it approximates the gamma function almost as well, (it needs to be checked).
example thinBasic script, look carefully for the changes. :)
Uses "console"
Uses "math"
Const n = 6
Dim As Ext g = 5
Dim As Ext a(n,n), b(n,n), c(n,1)
Dim As Integer i, j
'for n = 6, the following sets up matrix a as:
' [[ 1, 1/2, 1/3, 1/4, 1/5, 1/6 ]
' [ 1/2, 1/3, 1/4, 1/5, 1/6, 1/7 ]
' [ 1/3, 1/4, 1/5, 1/6, 1/7, 1/8 ]
' [ 1/4, 1/5, 1/6, 1/7, 1/8, 1/9 ]
' [ 1/5, 1/6, 1/7, 1/8, 1/9, 1/10 ]
' [ 1/6, 1/7, 1/8, 1/9, 1/10, 1/11 ]]
For i=1 To n
'<< note the change
For j=1 To n '<< note the change
a(i,j)=1/(i+j-1) '<< note the change
Next
Next
For i=1 To n
c(i,1) = Factorial(i-1)/((i+g-0.5)^(i-0.5)/Exp(i+g-0.5)*Sqr(2*M_PI))-1 '<< note the change
Next
MAT b() = INV(a())
MAT a() = b() * c()
Console_WriteLine "C( 1) = 1"
For i=1 To n
Console_WriteLine "C("&Str$(i+1)&") = "&a(i,1)
Next
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
here are the coefficients as calculated with Maple
c[0] := 1
c[1] := 76.180091728331374539146745
c[2] := -86.505320289513654614503042
c[3] := 24.014097921606006186204527
c[4] := -1.2317386147754424755938420
c[5] := .120745388047503859955038e-2
c[6] := -.4868518292851569517614e-5
danbaron
15-08-2011, 03:53
It seems that the new values disagree more with the values from Numerical Recipes.
' code -----------------------------------------------------------------------------
Uses "console"
Uses "math"
Const n = 6
Dim As Ext a(n), b(n), c(n)
Dim As Integer i
' From Johan.
a(1) = 76.1800917283313588
a(2) = -86.5053202895132074
a(3) = 24.0140979216034793
a(4) = -1.23173861476925595
a(5) = 1.2074538708475302 * 10^-3
a(6) = -4.8685153615224408 * 10^-6
' From Numerical Recipes.
b(1) = 76.18009172947146
b(2) = -86.50532032941677
b(3) = 24.01409824083091
b(4) = -1.231739572450155
b(5) = 0.1208650973866179 * 10^-2
b(6) = -0.5395239384953 * 10^-5
' abs(difference).
For i = 1 To n
c(i) = Abs(a(i)-b(i))
Next
For i = 1 To n
PrintL c(i)
Next
Console_WaitKey
' output ---------------------------------------------------------------------------
1.14010120094709677E-9
3.99035626000299182E-8
3.19227430699633996E-7
9.57680899050031287E-7
1.19710301864880011E-6
5.267240234305592E-7