Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]

  • To: mathgroup at smc.vnet.net
  • Subject: [mg84299] Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
  • From: Andrew Moylan <andrew.j.moylan at gmail.com>
  • Date: Tue, 18 Dec 2007 02:19:10 -0500 (EST)
  • References: <fk73ks$65f$1@smc.vnet.net>

Due to line-wrapping issues, it might not work to directly paste the
definitions for the matrix m and vector b into Mathematica. Sorry
about that. I've put the definitions of m and b into a Mathematica
notebook that you can get from http://andrew.j.moylan.googlepages.com/files.
This might be more convenient.

Andrew Moylan


On Dec 18, 11:18 am, "Andrew Moylan" <andrew.j.moy... at gmail.com>
wrote:
> At the bottom of this message I define a machine-precision matrix m and a
> machine precision vector b. (I apologise for their large size.)
>
> 1. LinearSolve[m, b] is not equivalent to LinearSolve[m][b].
>
> The documentation for LinearSolve states that "LinearSolve[m,b] is
> equivalent to LinearSolve[m][b]." For m and b, however:
>
> sol1 = LinearSolve[m, b]
>
> {-1.68129*10^9,-7.29642*10^11,-6.31285*10^20,4.99506*10^22,3.84561*10^24,-1.
> 3535*10^26,-7.82045*10^25,1.69903*10^27,2.25713*10^25,-2.49272*10^26,-1.0616
> *10^26,9.03507*10^26,-4.77786*10^24,2.87158*10^25,4.55334*10^25,-2.43647*10^
> 26,-2.82578*10^26,1.27101*10^27,-1.38448*10^28,5.52681*10^30}
>
> sol2 = LinearSolve[m][b]
>
> {-1.8893*10^-6,0.00627488,-0.000025753,-0.00767333,0.0000624517,-0.000027213
> 7,-0.0000464774,0.00188393,7.51337*10^-6,-0.0000585884,6.01809*10^-6,-0.0004
> 30135,-6.4536*10^-7,-0.0000946245,-1.29985*10^-6,0.0000701392,-4.03265*10^-7
> ,0.0000499027,4.97681*10^-7,9.8366*10^-6}
>
> Needless to say, sol2 is the one that approximately solves m.x==b. Type
> (m.sol1 - b) and (m.sol2 - b) to see this. sol1 is wild.
>
> 2. When using arbitrary precision arithmetic instead of machine precision
> arithmetic, the wild behaviour does not occur.
>
> Here are m and b at precision $MachinePrecision (i.e., 15.9546):
> mm = SetPrecision[m, $MachinePrecision];
> bb = SetPrecision[b, $MachinePrecision];
>
> For these definitions, LinearSolve[mm,bb] and LinearSolve[mm][bb] *are*
> approximately equal, and they are approximately equal to sol2 (the tame
> solution).
>
> Can anyone shed some light on this topic? My working hypothesis is that 1.
> sol1 is the result of a bug (or inconsistent/bad handling of poorly
> conditioned matrices); and 2. that bug isn't present in the arbitrary
> precision version of LinearSolve. Then the workaround would be: don't use
> LinearSolve[m, b]; instead use LinearSolve[m][b].
>
> Andrew Moylan
>
> Here are the definitions of m and b:
>
> m={{0.,120.,4.,-120.,-16.,120.,36.,-120.,-64.,120.,100.,-120.,-144.,120.,196
> .,-120.,-256.,120.,324.,-120.},{-120.,0.,120.,4.,-120.,-16.,120.,36.,-120.,-
> 64.,120.,100.,-120.,-144.,120.,196.,-120.,-256.,120.,324.},{0.,137.667770123
> 3155,4.,-126.52361302070159,-14.704805681953697,94.89537283269993,28.5433706
> 5200333,-47.90363918197348,-40.5445084603741,-6.84365962903201,45.5838278268
> 89575,60.48297480676686,-39.728436874545594,-104.33013932096895,21.378320417
> 681632,131.28632837898377,8.062099516267708,-136.987354117796,-44.1577316683
> 6576,120.51022451941246},{-137.6677701233155,0.,126.52361302070159,4.,-94.89
> 537283269993,-14.704805681953697,47.90363918197348,28.54337065200333,6.84365
> 962903201,-40.5445084603741,-60.48297480676686,45.583827826889575,104.330139
> 32096895,-39.728436874545594,-131.28632837898377,21.378320417681632,136.9873
> 54117796,8.062099516267708,-120.51022451941246,-44.15773166836576},{0.,190.5
> 9730297009568,4.,-140.79278017868444,-11.819078485136764,17.407811913964213,
> 14.191990544591762,115.0747432704792,-4.317884714201286,-187.4174897103089,-
> 15.679345960738772,161.81300927620237,34.27414014713948,-51.64259235864097,-
> 37.12435871470173,-85.51702774336465,16.98339066468509,177.98415019794658,19
> .503939521116326,-177.43405229412292},{-190.59730297009568,0.,140.7927801786
> 8444,4.,-17.407811913964213,-11.819078485136764,-115.0747432704792,14.191990
> 544591762,187.4174897103089,-4.317884714201286,-161.81300927620237,-15.67934
> 5960738772,51.64259235864097,34.27414014713948,85.51702774336465,-37.1243587
> 1470173,-177.98415019794658,16.98339066468509,177.43405229412292,19.50393952
> 1116326},{0.,322.872632447701,4.,-154.48510511279565,-7.655531727995272,-175
> .03942948598953,-1.0111563678108242,321.98759337281865,16.60122006528058,-13
> 3.08385014981712,-18.17273422029246,-194.634138709607,-4.033538586739716,319
> .3373281801208,29.944998291495985,-110.95299226741214,-27.37121776748681,-21
> 3.1618091027223,-9.033977175382377,314.93636636525736},{-322.872632447701,0.
> ,154.48510511279565,4.,175.03942948598953,-7.655531727995272,-321.9875933728
> 1865,-1.0111563678108242,133.08385014981712,16.60122006528058,194.6341387096
> 07,-18.17273422029246,-319.3373281801208,-4.033538586739716,110.952992267412
> 14,29.944998291495985,213.1618091027223,-27.37121776748681,-314.936366365257
> 36,-9.033977175382377},{0.,688.6670919150857,4.,-114.72119024973249,-2.66535
> 03057498296,-650.4455719765721,-10.667982702682368,331.4293532674049,10.0696
> 8577122105,540.0236559683539,13.586339913460915,-511.348380335829,-20.536381
> 943186196,-369.6583357342561,-11.038447478035566,634.5070001073631,31.584889
> 64242715,158.26040739218521,2.3538027804037114,-687.2344282612197},{-688.667
> 0919150857,0.,114.72119024973249,4.,650.4455719765721,-2.6653503057498296,-3
> 31.4293532674049,-10.667982702682368,-540.0236559683539,10.06968577122105,51
> 1.348380335829,13.586339913460915,369.6583357342561,-20.536381943186196,-634
> .5070001073631,-11.038447478035566,-158.26040739218521,31.58488964242715,687
> .2344282612197,2.3538027804037114},{0.,1955.742698127284,4.,319.116359538738
> 95,2.610702193856547,-1851.6029756343717,-10.722043885311114,-923.3643533687
> 346,-9.886748303381077,1550.2743077507857,13.837048847537616,1429.2774204092
> 866,20.248784522486616,-1083.847095614267,-11.662601235917451,-1782.97766919
> 99382,-31.348033453698886,501.9941315420778,3.4859696207965745,1946.79731676
> 494},{-1955.742698127284,0.,-319.11635953873895,4.,1851.6029756343717,2.6107
> 02193856547,923.3643533687346,-10.722043885311114,-1550.2743077507857,-9.886
> 748303381077,-1429.2774204092866,13.837048847537616,1083.847095614267,20.248
> 784522486616,1782.9776691999382,-11.662601235917451,-501.9941315420778,-31.3
> 48033453698886,-1946.79731676494,3.4859696207965745},{0.,8267.868082849189,4
> .,3927.69394068881,7.6008836161019895,-4536.125017493609,-1.1674815477135652
> ,-8237.513731446,-16.680749126890216,-3290.4228773517743,-17.86489003655355,
> 5111.248564132217,4.652781192989794,8146.673560968948,30.168278100627578,262
> 8.9911352802987,26.554230419975788,-5648.841605090274,-10.404652778807701,-7
> 996.014586041039},{-8267.868082849189,0.,-3927.69394068881,4.,4536.125017493
> 609,7.6008836161019895,8237.513731446,-1.1674815477135652,3290.4228773517743
> ,-16.680749126890216,-5111.248564132217,-17.86489003655355,-8146.67356096894
> 8,4.652781192989794,-2628.9911352802987,30.168278100627578,5648.841605090274
> ,26.554230419975788,7996.014586041039,-10.404652778807701},{0.,66714.5128785
> 9093,4.,49053.640102814905,11.764430373243481,5421.5038143728425,13.95034162
> 6292579,-41081.027085007,3.8241097110929165,-65833.36416473465,-16.221120559
> 06104,-55730.47653404341,-34.361001075598665,-16121.29969181952,-36.24174857
> 068484,32023.238065218116,-15.09440727983296,63213.194009826875,21.624753701
> 270276,60935.1643846494},{-66714.51287859093,0.,-49053.640102814905,4.,-5421
> .5038143728425,11.764430373243481,41081.027085007,13.950341626292579,65833.3
> 6416473465,3.8241097110929165,55730.47653404341,-16.22112055906104,16121.299
> 69181952,-34.361001075598665,-32023.238065218116,-36.24174857068484,-63213.1
> 94009826875,-15.09440727983296,-60935.1643846494,21.624753701270276},{0.,2.1
> 09193656453301*^6,4.,1.931251213300796*^6,14.650157570060415,1.4274481663250
> 76*^6,28.242584405174718,682791.3566687256,39.659403480135936,-177073.045991
> 05718,43.712855697405764,-1.0070598598161684*^6,36.57092835044224,-1.6671251
> 576076704*^6,16.935107017831488,-2.0458959211793535*^6,-13.318096903865762,-
> 2.0794620445450086*^6,-49.21136470378096,-1.7621599055136922*^6},{-2.1091936
> 56453301*^6,0.,-1.931251213300796*^6,4.,-1.427448166325076*^6,14.65015757006
> 0415,-682791.3566687256,28.242584405174718,177073.04599105718,39.65940348013
> 5936,1.0070598598161684*^6,43.712855697405764,1.6671251576076704*^6,36.57092
> 835044224,2.0458959211793535*^6,16.935107017831488,2.0794620445450086*^6,-13
> .318096903865762,1.7621599055136922*^6,-49.21136470378096},{0.,3.21113804748
> 0835*^10,4.,3.2001703830231155*^10,15.945351888106718,3.1673423099192696*^10
> ,35.67267128166534,3.1128780771958237*^10,62.91151228627201,3.03714973078828
> 5*^10,97.28714180413428,2.94067457209015*^10,138.32438810885077,2.8241116242
> 59851*^10,185.4525285933514,2.6882571304243427*^10,238.01117208726953,2.5340
> 39114531167*^10,295.2570820108493,2.362511042003884*^10},{-3.211138047480835
> *^10,0.,-3.2001703830231155*^10,4.,-3.1673423099192696*^10,15.94535188810671
> 8,-3.1128780771958237*^10,35.67267128166534,-3.037149730788285*^10,62.911512
> 28627201,-2.94067457209015*^10,97.28714180413428,-2.824111624259851*^10,138.
> 32438810885077,-2.6882571304243427*^10,185.4525285933514,-2.534039114531167*
> ^10,238.01117208726953,-2.362511042003884*^10,295.2570820108493}};
>
> b={1.4715177646857693,0.,1.7272858774111686,0.,2.4384494832359858,0.,3.86819
> 71196755,0.,6.1130185623361015,0.,7.4508820294899945,0.,3.393321220405215,0.
> ,0.03385691675272283,0.,3.3743307736651725*^-14,0.,1.23885112322841128107475
> 3`10.061735978331926*^-496,0.};



  • Prev by Date: Re: Table[Plot[]]doesn't work as it should in v6
  • Next by Date: Re: Updating NumberPadding?
  • Previous by thread: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
  • Next by thread: Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]