Results 1 to 10 of 12

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

Hybrid View

Previous Post Previous Post   Next Post Next Post
  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.

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
  •