# An Interesting, and Perhaps New, Matrix

I do not recall having seen this matrix before, but I will not be surprised to learn that somebody else already knows all about it, especially if that person's name is Nick.

### Contents

#### Q

I've been investigating the matrices generated by this elegant one-liner.

    Q = @(n) (-n:n).^2 + (-n:n)'.^2;


The Q is for "quadratic".

The middle column contains the squares of the integers from -n to n. So does the middle row. The apostrophe summons singleton expansion. The resulting matrix has order 2*n+1. Here is Q(5).

    Q5 = Q(5)

Q5 =

50    41    34    29    26    25    26    29    34    41    50
41    32    25    20    17    16    17    20    25    32    41
34    25    18    13    10     9    10    13    18    25    34
29    20    13     8     5     4     5     8    13    20    29
26    17    10     5     2     1     2     5    10    17    26
25    16     9     4     1     0     1     4     9    16    25
26    17    10     5     2     1     2     5    10    17    26
29    20    13     8     5     4     5     8    13    20    29
34    25    18    13    10     9    10    13    18    25    34
41    32    25    20    17    16    17    20    25    32    41
50    41    34    29    26    25    26    29    34    41    50



I like the contour plot.

    contourf(Q(100))
axis square
colorbar
title('Q(100)')


#### D

For another blog post under development, I need a logical mask that carves a circular region out of graphic. This disc does the job.

    D = @(n) Q(n) <= n^2;


Here is my carver.

    spy(D(100))
title('D(100)')


Did you notice the digits in the count of nonzeros in D(100)? It happens whenever n is a power of 10.

    fprintf('%15s %12s\n','n','nnz(D(n))')
for n = 10.^(0:4)
fprintf('%15d %12d\n',n, nnz(D(n)))
end

              n    nnz(D(n))
1            5
10          317
100        31417
1000      3141549
10000    314159053


#### O.E.I.S.

A classic problem, described in the Online Encyclopedia of Integer Sequences, asks how many points with integer coordinates lie within a disc of increasing radius. Our nonzero count provides the answer.

    fprintf('%15s %8s\n','n','a(n)')
for n = [1:15  99:101  499:501  999:1001]
if mod(n,100) == 99
fprintf('%15s %8s\n','-','-')
end
a(n) = nnz(D(n));
fprintf('%15d %8d\n',n,a(n))
end

              n     a(n)
1        5
2       13
3       29
4       49
5       81
6      113
7      149
8      197
9      253
10      317
11      377
12      441
13      529
14      613
15      709
-        -
99    30757
100    31417
101    32017
-        -
499   782197
500   785349
501   788509
-        -
999  3135157
1000  3141549
1001  3147833


#### R

Taking the reciprocals of the matrix entries, and reducing the range of the anonymous index, produces a matrix that behaves a bit like the Hilbert matrix, hilb(n).

    R = @(n) 1./((1:n).^2 + (1:n)'.^2);


Here are the 5-by-5's.

    format rat
R5 = R(5)
H5 = hilb(5)

R5 =

1/2            1/5            1/10           1/17           1/26
1/5            1/8            1/13           1/20           1/29
1/10           1/13           1/18           1/25           1/34
1/17           1/20           1/25           1/32           1/41
1/26           1/29           1/34           1/41           1/50

H5 =

1              1/2            1/3            1/4            1/5
1/2            1/3            1/4            1/5            1/6
1/3            1/4            1/5            1/6            1/7
1/4            1/5            1/6            1/7            1/8
1/5            1/6            1/7            1/8            1/9



#### Condition

Going away from the diagonal, the elements of R(n) decay more rapidly than those of hilb(n), so R(n) is better conditioned than hilb(n).

    format short e
fprintf('%15s %12s %12s\n','n','cond R','cond H')
for n = 1:12
fprintf('%15d %12.2e %12.2e\n',n,cond(R(n)),cond(hilb(n)))
end

              n       cond R       cond H
1     1.00e+00     1.00e+00
2     1.53e+01     1.93e+01
3     2.04e+02     5.24e+02
4     2.59e+03     1.55e+04
5     3.21e+04     4.77e+05
6     3.89e+05     1.50e+07
7     4.67e+06     4.75e+08
8     5.54e+07     1.53e+10
9     6.53e+08     4.93e+11
10     7.65e+09     1.60e+13
11     8.92e+10     5.22e+14
12     1.04e+12     1.62e+16


#### Extra Credit

What is the rank of Q(n)? Why? See a paper in SIAM Review by Strang and Moler.

Why is the table of values for nnz(D(10^k)) so short? How might you extend this table?

Investigate R(n). Is it positive definite? What are its eigenvalues? What is its inverse? What is the sign pattern of the elements of its inverse? For what values of n can you compute the inverse reliably using floating point arithmetic? How does all this compare with hilb(n) and invhilb(n) ?

Comments in the Comments, or in email to me, are welcome.

Published with MATLAB® R2022a

|

### 댓글

댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.