exponenta event banner

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

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

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

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

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

H = hilb(20);

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

cond(H)
ans = 2.1065e+18

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

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

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

    1.0000   -0.0000   -0.0003    0.0044    0.0212   -0.2436    0.4383   -0.6683   -0.4409    0.6464   -0.9986   -0.4477    8.3067   -0.2302   -0.4680   -0.4770    1.4327   -1.4875    8.0000    2.0000
    0.4314    1.0000   -0.0002    0.0003    0.0528   -0.0284    0.9397    0.6603   -0.1158    1.3085   -0.5132    0.2835    0.5888   -1.7213    0.7195   -1.7957    0.6355    2.8214   -8.0000    2.0000
    0.7987   -0.8048    0.9997   -0.0041    0.1027    0.1222    0.5629   -0.1333    0.5837    1.0126   -0.9596    0.1601   -0.3606   -2.5592    0.5373   -2.9668   -0.2067    3.8693   -8.0000         0
    1.0486   -1.6648    0.4041    0.9885    0.0958   -0.1427    1.3149    0.2797   -0.1724    0.8018   -0.1066    0.3468    5.6058   -2.5304    2.5296   -3.6839   -0.4191    3.5598  -16.0000    1.0000
    1.2164   -2.4406    1.0946   -0.1853    1.1204   -0.1147    1.2598   -0.9847    0.2424   -0.0230    0.5197    0.8731    3.8277   -4.7249    0.5983   -2.3525   -0.3189    3.4062   -4.0000    3.0000
    1.3286   -3.1054    1.9447   -0.5911    0.1962    0.8550    1.2190   -1.0755    0.2927    1.1158    0.2254    0.9657    3.0420   -2.0190    1.2346   -1.7681   -0.4577    3.5038         0         0
    1.4027   -3.6617    2.8529   -1.1814    0.3594   -0.2091    2.4910   -0.4949    0.6416    0.4732    0.2943   -0.5944    4.9283   -4.1578    1.6593   -3.0688   -0.1768    4.6058         0    2.0000
    1.4501   -4.1206    3.7522   -1.8860    0.5367    0.0995    0.2101    1.4302   -1.0067    1.5837   -1.5936    3.3039   -0.4865    0.0223    0.9902   -1.7856   -0.3147    2.2500         0         0
    1.4784   -4.4954    4.6045   -2.6478    0.8130    0.1115    0.0798    0.1202    1.1006    0.0361   -0.2599    3.4626   -0.4334   -0.7285   -0.3979    0.2157    0.6016    0.2411         0   -1.0000
    1.4930   -4.7986    5.3893   -3.4322    1.1685   -0.0259    0.4050   -0.0399    1.2269    0.7081    0.0601    0.9479    3.7747   -2.1672    1.9212   -1.5514    0.2045    2.7134   -4.0000         0
      ⋮

Теперь используйте 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.0004    0.0013   -0.0037    0.0047   -0.2308    0.4019    0.2620   -0.4443   -7.5099    3.4753    3.2884   -1.1618    0.2301    0.4295    0.0537
   -0.0000    0.0000   -0.0000    0.0000   -0.0000   -0.0001   -0.0009    0.0172   -0.0628    0.1251   -0.8975    2.7711   -2.8924   -1.4888   -2.5102    6.4897   -3.0581    0.7925         0   -0.0268
    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0001    0.0009    0.0042   -0.0303    0.0271   -0.1028    0.2862   -0.9267   -4.0597    1.8194    4.7727   -1.3255    0.4667    0.2147   -0.0268
    0.0000    0.0000   -0.0000   -0.0000    0.0000   -0.0001    0.0004    0.0002   -0.0056    0.0526   -0.1136    0.3616   -1.6920   -3.7214    0.2709    4.4169   -1.1602    0.5541         0   -0.0268
   -0.0000    0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0006    0.0068   -0.0394    0.1449   -0.7818    1.9314   -2.1273   -1.0745   -2.2988    4.6349   -1.8923    0.3793    0.2147    0.0537
    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0014   -0.0028    0.0304   -0.1053   -0.0724   -0.2073   -0.3769   -5.0729    2.7552    2.3488   -0.5046    0.2240         0         0
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0013    0.0101   -0.0520    0.1307   -0.7715    2.9297   -3.3374    1.4084   -3.8703    7.1053   -2.7578    0.8815   -0.2147         0
    0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0002    0.0018   -0.0040    0.0231   -0.0892    0.2000    0.2766   -0.8262   -4.7229    1.6513    2.1857   -0.4749   -0.1747         0         0
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000    0.0004    0.0144   -0.0513    0.1165   -0.7063    1.9776   -2.0602   -1.0763   -3.2081    3.8539   -2.4580    0.5675         0    0.0268
    0.0000    0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0008   -0.0019    0.0027   -0.0167   -0.0497    0.8865   -0.6720   -1.2729    0.0571    2.3625   -1.2616    0.2443    0.2147         0
      ⋮

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

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)]