Гильбертовы матрицы и их инверсии

В этом примере показано, как вычислить инверсию Гильбертовой матрицы с помощью Symbolic Math Toolbox™.

Определение: Гильбертова матрица является квадратной матрицей с записями, являющимися модульной частью. Hij=1i+j-1. Например, 3x3 Гильбертова матрица H=[11213121314131415]

Символьные расчеты дают точные результаты для этих плохо обусловленных матриц, в то время как чисто численные методы перестали работать.

Создайте 20 20 числовую Гильбертову матрицу.

H = hilb(20);

Найдите число обусловленности этой матрицы. Гильбертовы матрицы плохо обусловлены, означая, что у них есть большие числа обусловленности, указывающие, что такие матрицы почти сингулярны. Обратите внимание на то, что вычислительные числа обусловленности также подвержены числовым ошибкам.

cond(H)
ans = 2.1211e+18

Поэтому инвертирование Гильбертовых матриц численно неустойчиво. Когда вы вычисляете обратную матрицу, H*inv(H) должен возвратить единичную матрицу или матрицу близко к единичной матрице в некотором допуске на погрешность.

Во-первых, вычислите инверсию H при помощи inv функция. Предупреждение выдано из-за числовой нестабильности.

H*inv(H)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.549910e-19.
ans = 20×20

    1.0000    0.0000   -0.0000   -0.0024   -0.0156    0.0625    0.0312   -0.2500    0.8320         0    0.0312    1.3750    0.1250   -0.3438   -0.8750   -0.5000         0         0   -1.2500         0
   -0.1868    1.0000   -0.0001   -0.0024    0.0273    0.0781   -0.5625   -0.5625   -1.0273    0.3125   -0.4082   -0.6250    3.3750   -0.4219    0.5000   -0.1250   -2.0000    1.5000   -1.7500   -0.1250
   -0.2244    0.0051    1.0001   -0.0054    0.0273    0.0625   -0.1562    0.8750   -1.0781    0.2500   -0.1113    1.8750    1.7500    0.6250    0.1250   -2.8750   -2.0000   -0.5000    0.7500    0.6875
   -0.2236   -0.0412    0.0613    0.9985    0.0273    0.0625    0.3438   -0.2500    0.2637    0.3750   -0.8867    1.2500   -1.5000   -0.3906   -0.8750    0.5000   -3.0000    1.0000    3.5000   -0.0625
   -0.2112   -0.1026    0.1511   -0.0222    1.0273    0.1250    0.0938    0.1250    0.7441   -0.5625   -0.2285    0.1250    0.7500   -1.7188    0.3125   -2.0000   -3.5000   -2.0000    1.5000   -0.5000
   -0.1962   -0.1541    0.1993    0.0381   -0.0254    0.9531    0.0625    0.2500   -1.0879   -0.8125   -0.1777   -0.6250   -0.6250   -0.4688   -0.6875   -3.5000    0.2500    0.5000         0    0.7500
   -0.1815   -0.1888    0.1856    0.2017   -0.1602    0.0781    0.6250    1.3125    0.5137    0.8750    0.4668    2.2500   -2.2500   -0.6562    0.9375   -2.7500         0    2.0000    1.5000   -0.3125
   -0.1681   -0.2076    0.1196    0.4419   -0.3184   -0.2344    0.7812    0.4375    0.5254    0.5625   -0.0957   -0.6250   -0.1250   -1.0000   -1.4375    2.0000   -5.5000    2.7500    0.2500    0.4688
   -0.1563   -0.2138    0.0152    0.7461   -0.5840    0.0312   -0.2500    1.3438   -0.2051    2.3125    0.0449    1.3750   -8.0000    2.4922    2.8750   -5.0000    3.0000   -1.2500    1.7500    0.0938
   -0.1460   -0.2108   -0.1120    1.0696   -0.7773    0.0781    0.2500   -0.1250    0.8945    0.3125   -0.4688    1.5000   -5.0000   -1.1953    2.0625   -2.5625    2.2500   -1.7500    0.5000    0.1250
      ⋮

Теперь используйте MATLAB® invhilb функция, которая предлагает лучшую точность для Гильбертовых матриц. Эта функция находит точные инверсии Гильбертовых матриц до 15 15. Для 20 20 Гильбертовой матрицы, invhilb находит приближение обратной матрицы.

H*invhilb(20)
ans = 20×20
1010 ×

    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0003    0.0027   -0.0113   -0.0526   -0.1057    0.3087    0.3288   -0.3020   -7.1001    1.5435    4.8989   -0.1242   -0.0319    0.1770    0.0029
   -0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0001   -0.0015    0.0114   -0.0749    0.2441   -0.2819    3.6574   -4.3487   -1.2214   -4.3218    7.4826   -2.6407    0.8623   -0.0776    0.0029
    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0003   -0.0005   -0.0029    0.0180   -0.1460   -0.5184    1.0100    0.2483   -6.2948   -4.8050    5.7982   -1.0939    0.3758    0.0495   -0.0033
    0.0000   -0.0000    0.0000   -0.0000   -0.0000    0.0000    0.0012    0.0152   -0.0434   -0.1653   -0.2097    2.0099   -2.7850   -3.4360    0.9328    5.2613   -1.8220    0.6610   -0.0126   -0.0050
   -0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0001    0.0007    0.0127   -0.1088    0.1946   -0.4312    1.4630   -2.5032   -2.9998    1.2885    6.0666   -1.9898    0.3238   -0.0839    0.0045
    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0003   -0.0044    0.0231    0.0562         0    0.5838   -0.5302   -3.7447    4.9325    2.3287   -1.3220   -0.1242    0.1323   -0.0138
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0001   -0.0003    0.0170   -0.0938    0.0721   -0.5335    3.4326   -5.1540   -1.0536    0.5570    7.6974   -1.5234    0.2382   -0.1527   -0.0038
    0.0000    0.0000    0.0000   -0.0000    0.0000   -0.0003    0.0024   -0.0078    0.0091   -0.2278    0.5872    0.3792    0.7046   -4.3151    3.8319   -0.5369    0.1745   -0.0042    0.1141   -0.0141
   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0006    0.0113   -0.0708    0.1850   -0.7147    2.1643   -1.7851   -3.5970   -0.0201    4.1607   -1.3925    0.5637   -0.0497    0.0036
    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0001    0.0001    0.0023   -0.0149   -0.0679    0.3087    0.1074   -0.1745   -4.3218    0.8925    3.9259   -0.7986    0.2533    0.0455   -0.0091
      ⋮

Чтобы избежать ошибок округления, используйте точные символьные расчеты. Для этого создайте символьную Гильбертову матрицу.

Hsym = sym(H)
Hsym = 

(11213141516171819110111112113114115116117118119120121314151617181911011111211311411511611711811912012113141516171819110111112113114115116117118119120121122141516171819110111112113114115116117118119120121122123151617181911011111211311411511611711811912012112212312416171819110111112113114115116117118119120121122123124125171819110111112113114115116117118119120121122123124125126181911011111211311411511611711811912012112212312412512612719110111112113114115116117118119120121122123124125126127128110111112113114115116117118119120121122123124125126127128129111112113114115116117118119120121122123124125126127128129130112113114115116117118119120121122123124125126127128129130131113114115116117118119120121122123124125126127128129130131132114115116117118119120121122123124125126127128129130131132133115116117118119120121122123124125126127128129130131132133134116117118119120121122123124125126127128129130131132133134135117118119120121122123124125126127128129130131132133134135136118119120121122123124125126127128129130131132133134135136137119120121122123124125126127128129130131132133134135136137138120121122123124125126127128129130131132133134135136137138139)[sym (1), sym (1/2), sym (1/3), sym (1/4), sym (1/5), sym (1/6), sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20); sym (1/2), sym (1/3), sym (1/4), sym (1/5), sym (1/6), sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21); sym (1/3), sym (1/4), sym (1/5), sym (1/6), sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22); sym (1/4), sym (1/5), sym (1/6), sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23); sym (1/5), sym (1/6), sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24); sym (1/6), sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25); sym (1/7), sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26); sym (1/8), sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27); sym (1/9), sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28); sym (1/10), sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29); sym (1/11), sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30); sym (1/12), sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31); sym (1/13), sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32); sym (1/14), sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33); sym (1/15), sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33), sym (1/34); sym (1/16), sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33), sym (1/34), sym (1/35); sym (1/17), sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33), sym (1/34), sym (1/35), sym (1/36); sym (1/18), sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33), sym (1/34), sym (1/35), sym (1/36), sym (1/37); sym (1/19), sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33), sym (1/34), sym (1/35), sym (1/36), sym (1/37), sym (1/38); sym (1/20), sym (1/21), sym (1/22), sym (1/23), sym (1/24), sym (1/25), sym (1/26), sym (1/27), sym (1/28), sym (1/29), sym (1/30), sym (1/31), sym (1/32), sym (1/33), sym (1/34), sym (1/35), sym (1/36), sym (1/37), sym (1/38), sym (1/39)]

Получите значение числа обусловленности. Это было выведено символьными методами и свободно от числовых ошибок.

vpa(cond(Hsym))
ans = 24521565858153031724608315432.509vpa ('24521565858153031724608315432.509')

Несмотря на то, что его число обусловленности является большим, можно вычислить точную инверсию матрицы.

Hsym*inv(Hsym)
ans = 

(1000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001)[sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1), sym (0); sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (0), sym (1)]

Для просмотра документации необходимо авторизоваться на сайте