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

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

Определение: Гильбертова матрица является квадратной матрицей с записями, являющимися модульной частью. $H_{i j} = \frac{1}{i + j - 1}$ . Например, Гильбертова матрица 3x3 $H = [\begin{array}{ccc} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{array}]$

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

Создайте 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
10¹⁰ ×

    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 = 
 $(\begin{matrix} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} \\ \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} \\ \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} \\ \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} \\ \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} \\ \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} \\ \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} \\ \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} \\ \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} \\ \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} \\ \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} \\ \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} \\ \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} \\ \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} & \frac{1}{37} \\ \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} & \frac{1}{37} & \frac{1}{38} \\ \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} & \frac{1}{37} & \frac{1}{38} & \frac{1}{39} \end{matrix})$

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

vpa(cond(Hsym))

ans =  $24521565858153031724608315432.509$

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

Hsym*inv(Hsym)

ans = 
 $(\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix})$

Документация

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

Документация Symbolic Math Toolbox

Поддержка