Page 1 of 2 12 LastLast
Results 1 to 10 of 12

Thread: how do you solve a system of linear equations in ThinBasic?

  1. #1

    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.
    Last edited by jack; 14-08-2011 at 01:27.

  2. #2
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    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
    
    Last edited by danbaron; 13-08-2011 at 07:32.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  3. #3
    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
    
    Last edited by jack; 13-08-2011 at 12:58.

  4. #4
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    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
    
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  5. #5
    Super Moderator Petr Schreiber's Avatar
    Join Date
    Aug 2005
    Location
    Brno - Czech Republic
    Posts
    7,128
    Rep Power
    732
    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
    Learn 3D graphics with ThinBASIC, learn TBGL!
    Windows 10 64bit - Intel Core i5-3350P @ 3.1GHz - 16 GB RAM - NVIDIA GeForce GTX 1050 Ti 4GB

  6. #6
    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.


    [edit]
    fixed my first post.
    [/edit]
    Last edited by jack; 14-08-2011 at 01:29.

  7. #7
    Member REDEBOLT's Avatar
    Join Date
    Dec 2006
    Location
    Pickerington, Ohio, USA
    Age
    87
    Posts
    62
    Rep Power
    24
    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
    Attached Files Attached Files
    Last edited by REDEBOLT; 14-08-2011 at 07:06. Reason: posted the zip file.

  8. #8
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    (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.
    Last edited by danbaron; 14-08-2011 at 08:23.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  9. #9
    Member REDEBOLT's Avatar
    Join Date
    Dec 2006
    Location
    Pickerington, Ohio, USA
    Age
    87
    Posts
    62
    Rep Power
    24
    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 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

  10. #10
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    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.


    Last edited by danbaron; 14-08-2011 at 12:18.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

Page 1 of 2 12 LastLast

Similar Threads

  1. operating system in thinbasic architecture ?
    By ioptic in forum thinBasic General
    Replies: 8
    Last Post: 14-01-2009, 23:53
  2. New advertising system in thinBasic community forum
    By ErosOlmi in forum Web and Forum
    Replies: 3
    Last Post: 04-08-2008, 18:07

Members who have read this thread: 0

There are no members to list at the moment.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •