OK, thank you. I still think that there is something wrong with %.
I have tried the same sort of example as I first introduced here
with real numbers rather than complex numbers. And J still is wrong
on it by several orders of magnitude. Here is the program in J (which does
poorly, absolute error about 0.003) and in Octave (which does the same
essentially perfectly, absolute error 3e-15).
Both J and Octave do the QR decomposition, how can there be such
a wide difference?
Thank you very much!
Yours sincerely,
Imre Patyi
x=:_0.5+(i.1000)%999
b=:^x
X=:(1 o. (x */ >:i.10)) ,. (2 o. (x */ i.11))
,.c=:b %. X
0j15":>./|(X (+/ . *) c)-b
y=:_0.5+(i.1e6)%999999
0j15":>./|(((1 o. (y */ (>:i.10))) (+/ . *) ((i.10){c)) + ((2 o. (y */
(i.11))) (+/ . *) ((10+i.11){c)))-(,.^y)
NB. Sample run. Windows 10, i5, J805.
NB. x=:_0.5+(i.1000)%999
NB. b=:^x
NB. X=:(1 o. (x */ >:i.10)) ,. (2 o. (x */ i.11))
NB. ,.c=:b %. X
NB. 9.98124e8
NB. _1.40794e9
NB. 1.15836e9
NB. _6.34655e8
NB. 2.23796e8
NB. _3.80026e7
NB. _5.72319e6
NB. 4.98532e6
NB. _1.17017e6
NB. 105348
NB. 97984.1
NB. 1.97192e7
NB. _5.38456e7
NB. 6.1024e7
NB. _3.64392e7
NB. 7.47621e6
NB. 5.76293e6
NB. _5.64701e6
NB. 2.2945e6
NB. _487483
NB. 44448
NB. 0j15":>./|(X (+/ . *) c)-b
NB. 0.003757055930316
NB. y=:_0.5+(i.1e6)%999999
NB. 0j15":>./|(((1 o. (y */ (>:i.10))) (+/ . *) ((i.10){c)) + ((2 o. (y
*/ (i.11))) (+/ . *) ((10+i.11){c)))-(,.^y)
NB. 0.003757055988523
x=linspace(-0.5,0.5,1000)' ;
b=exp(x) ;
X=[sin(x.*(1:10)) , cos(x.*(0:10))] ;
c=X\b
norm(X*c-b,"inf")
y=linspace(-0.5,0.5,1e6)' ;
norm(sin(y.*(1:10))*c(1:10)+cos(y.*(0:10))*c(11:end)-exp(y),"inf")
%{
Sample run.
GNU Octave, version 5.2.0
Copyright (C) 2020 John W. Eaton and others.
Octave was configured for "x86_64-w64-mingw32".
>> x=linspace(-0.5,0.5,1000)' ;
>> b=exp(x) ;
>> X=[sin(x.*(1:10)) , cos(x.*(0:10))] ;
>> c=X\b
c =
1.915512522
-0.121069735
-0.756548522
0.760884087
-0.454428637
0.190450190
-0.057022210
0.011742016
-0.001500367
0.000090086
1.899201318
-0.404596102
-1.085720223
0.971184293
-0.559477240
0.243764657
-0.081572358
0.020470129
-0.003643514
0.000411202
-0.000022161
>> norm(X*c-b,"inf")
ans = 2.7756e-15
>> y=linspace(-0.5,0.5,1e6)' ;
>> norm(sin(y.*(1:10))*c(1:10)+cos(y.*(0:10))*c(11:end)-exp(y),"inf")
ans = 2.9976e-15
>>
%}
On Thu, Apr 29, 2021 at 2:29 PM Henry Rich <[email protected]> wrote:
> %. does not use SVD. It does do QR decomposition followed by inverse.
>
> You can call LAPACK to perform the SVD, but you might need to delete
> columns with low singular values.
>
> Which routine you need to use depends on details of what you are going
> to do.
>
> Henry Rich
>
>
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
