% sperimentazione numerica su sistemi lineari rettangolari % sistema sovradeterminato a rango pieno A*x=b A=rand(10,7) A = 0.3517 0.0759 0.1622 0.4505 0.1067 0.4314 0.8530 0.8308 0.0540 0.7943 0.0838 0.9619 0.9106 0.6221 0.5853 0.5308 0.3112 0.2290 0.0046 0.1818 0.3510 0.5497 0.7792 0.5285 0.9133 0.7749 0.2638 0.5132 0.9172 0.9340 0.1656 0.1524 0.8173 0.1455 0.4018 0.2858 0.1299 0.6020 0.8258 0.8687 0.1361 0.0760 0.7572 0.5688 0.2630 0.5383 0.0844 0.8693 0.2399 0.7537 0.4694 0.6541 0.9961 0.3998 0.5797 0.1233 0.3804 0.0119 0.6892 0.0782 0.2599 0.5499 0.1839 0.5678 0.3371 0.7482 0.4427 0.8001 0.1450 0.2400 % termine noto b=rand(10,1) b = 0.4173 0.0497 0.9027 0.9448 0.4909 0.4893 0.3377 0.9001 0.3692 0.1112 % matrice A^T*A delle equazioni normali G=A'*A G = 3.9902 2.6847 2.8928 2.6795 3.1817 2.7188 2.1754 2.6847 2.4444 1.5821 2.0454 2.0509 1.3595 1.3471 2.8928 1.5821 2.9549 2.3942 2.9120 2.1905 1.5752 2.6795 2.0454 2.3942 3.2860 2.4976 1.8404 1.4822 3.1817 2.0509 2.9120 2.4976 3.8343 1.9284 1.7924 2.7188 1.3595 2.1905 1.8404 1.9284 2.5728 1.6185 2.1754 1.3471 1.5752 1.4822 1.7924 1.6185 1.8326 % termine noto delle equazioni normali c=A'*b c = 2.9635 2.4281 2.2784 2.8971 2.2280 1.8111 1.7096 % soluzione delle equazioni normali col metodo di Gauss con pivoting x1=G\c x1 = 0.0404 0.6275 0.7066 0.4365 -0.5977 -0.3741 0.3783 % norma del residuo norm(A*x1-b) ans = 0.5311 % fattorizzazione QR [Q R]=qr(A) Q = -0.1760 -0.2012 0.2289 0.4782 -0.0187 0.1193 0.7323 -0.0232 0.2653 0.1603 -0.4159 -0.6323 0.1195 -0.2100 -0.3200 -0.2123 0.0462 0.1298 -0.4046 -0.1968 -0.2930 0.1715 0.0433 -0.1991 0.4882 0.4023 0.2378 -0.3134 -0.3070 -0.4405 -0.2752 0.5124 -0.4511 0.0630 -0.1023 -0.4699 0.3462 0.2343 -0.2092 -0.0666 -0.4592 0.3967 0.3950 -0.3183 -0.4084 0.1673 -0.0895 0.0776 0.4048 -0.0430 -0.1431 -0.0781 -0.4455 0.3578 -0.4112 0.1268 -0.1909 -0.5252 0.1855 -0.3394 -0.3791 0.0743 0.3127 0.2549 0.2980 -0.4827 -0.3018 -0.4406 -0.0600 0.2783 -0.3773 -0.0472 -0.1068 0.4599 0.2533 0.2868 -0.3830 0.5776 0.0179 -0.0734 -0.1905 -0.3056 -0.3400 -0.3598 0.3956 -0.2468 0.0520 0.0278 0.6307 -0.0814 -0.2843 -0.0562 -0.3856 -0.2200 -0.0563 0.3714 0.0151 -0.1384 -0.1669 0.7284 R = -1.9976 -1.3440 -1.4482 -1.3414 -1.5928 -1.3611 -1.0890 0 0.7988 -0.4560 0.3037 -0.1124 -0.5881 -0.1460 0 0 -0.8061 -0.7320 -0.6874 0.0604 0.0849 0 0 0 0.9266 -0.1166 0.2563 0.1379 0 0 0 0 -0.8936 0.2621 -0.1296 0 0 0 0 0 -0.4863 -0.0901 0 0 0 0 0 0 0.7577 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 % soluzione delle equazioni normali mediante fattorizzazione QR R1=R(1:7,:) R1 = -1.9976 -1.3440 -1.4482 -1.3414 -1.5928 -1.3611 -1.0890 0 0.7988 -0.4560 0.3037 -0.1124 -0.5881 -0.1460 0 0 -0.8061 -0.7320 -0.6874 0.0604 0.0849 0 0 0 0.9266 -0.1166 0.2563 0.1379 0 0 0 0 -0.8936 0.2621 -0.1296 0 0 0 0 0 -0.4863 -0.0901 0 0 0 0 0 0 0.7577 c=Q'*b c = -1.4836 0.5436 -0.4688 0.4304 0.3871 0.1478 0.2866 0.0825 0.1154 -0.5118 c1=c(1:7) c1 = -1.4836 0.5436 -0.4688 0.4304 0.3871 0.1478 0.2866 c2=c(8:end) c2 = 0.0825 0.1154 -0.5118 x2=R1\c1 x2 = 0.0404 0.6275 0.7066 0.4365 -0.5977 -0.3741 0.3783 % confronto tra le soluzioni [x1 x2] ans = 0.0404 0.0404 0.6275 0.6275 0.7066 0.7066 0.4365 0.4365 -0.5977 -0.5977 -0.3741 -0.3741 0.3783 0.3783 x1-x2 ans = 1.0e-14 * 0.3920 -0.4219 -0.5218 0.1221 0.2887 0.0111 -0.0999 % soluzione mediante backslash (usa la fattorizzazione QR) x3=A\b x3 = 0.0404 0.6275 0.7066 0.4365 -0.5977 -0.3741 0.3783 [x1 x2 x3] ans = 0.0404 0.0404 0.0404 0.6275 0.6275 0.6275 0.7066 0.7066 0.7066 0.4365 0.4365 0.4365 -0.5977 -0.5977 -0.5977 -0.3741 -0.3741 -0.3741 0.3783 0.3783 0.3783 % soluzione mediante la matrice pseudoinversa Ap=pinv(A) Ap = 0.6070 0.1733 0.8118 -2.4886 0.8460 0.3329 -0.4454 1.4770 -1.1482 0.5726 -0.8380 -0.4302 -0.0235 1.8369 -0.0241 -0.6432 0.5531 -0.9663 0.8825 -0.2893 -0.4229 -0.2474 0.7516 0.5985 -0.6785 -0.3479 -0.4217 -0.2576 1.2495 0.6081 0.4589 -0.2937 -0.1238 -0.2053 -0.1891 0.5374 0.0441 0.6640 -0.5741 -0.0489 -0.2437 0.4740 -0.8515 0.3068 0.3796 0.4339 0.0371 -0.3556 -0.3075 -0.1651 -0.4245 0.4252 -0.8854 0.8815 -0.3221 -0.2139 1.0665 -0.4960 0.4948 -0.7674 0.9665 0.0610 0.3139 0.4570 -0.1181 -0.2520 -0.3983 -0.5055 0.0686 0.0199 x4=Ap*b x4 = 0.0404 0.6275 0.7066 0.4365 -0.5977 -0.3741 0.3783 [x1 x2 x3 x4] ans = 0.0404 0.0404 0.0404 0.0404 0.6275 0.6275 0.6275 0.6275 0.7066 0.7066 0.7066 0.7066 0.4365 0.4365 0.4365 0.4365 -0.5977 -0.5977 -0.5977 -0.5977 -0.3741 -0.3741 -0.3741 -0.3741 0.3783 0.3783 0.3783 0.3783 % sistema sottodeterminato a rango pieno A'*y=c A' ans = 0.3517 0.8308 0.5853 0.5497 0.9172 0.2858 0.7572 0.7537 0.3804 0.5678 0.0759 0.0540 0.5308 0.7792 0.9340 0.1299 0.5688 0.4694 0.0119 0.3371 0.1622 0.7943 0.3112 0.5285 0.1656 0.6020 0.2630 0.6541 0.6892 0.7482 0.4505 0.0838 0.2290 0.9133 0.1524 0.8258 0.5383 0.9961 0.0782 0.4427 0.1067 0.9619 0.0046 0.7749 0.8173 0.8687 0.0844 0.3998 0.2599 0.8001 0.4314 0.9106 0.1818 0.2638 0.1455 0.1361 0.8693 0.5797 0.5499 0.1450 0.8530 0.6221 0.3510 0.5132 0.4018 0.0760 0.2399 0.1233 0.1839 0.2400 c=rand(7,1) c = 0.7803 0.3897 0.2417 0.4039 0.0965 0.1320 0.9421 % soluzione delle equazioni normali del secondo tipo z=G\c z = 4.4552 -3.0021 -0.4896 0.8689 -1.2171 -2.9675 0.9614 y1=A*z y1 = 1.0611 -0.0516 0.8527 -0.5877 0.2933 -0.0817 -0.4469 0.4059 -0.3815 0.3628 % la soluzione verifica il sistema A'*y1-c ans = 1.0e-14 * -0.0999 0.0333 -0.0860 -0.0444 0.1721 0.0583 -0.0444 % soluzione mediante la matrice pseudoinversa y2=pinv(A')*c y2 = 1.0611 -0.0516 0.8527 -0.5877 0.2933 -0.0817 -0.4469 0.4059 -0.3815 0.3628 [y1 y2] ans = 1.0611 1.0611 -0.0516 -0.0516 0.8527 0.8527 -0.5877 -0.5877 0.2933 0.2933 -0.0817 -0.0817 -0.4469 -0.4469 0.4059 0.4059 -0.3815 -0.3815 0.3628 0.3628 norm(y1-y2) ans = 5.2093e-15 % soluzione mediante backslash (non e' di minima norma!) y3=A'\c y3 = 1.1498 -0.1788 0.7480 -0.6233 0.5500 0 -0.6800 0.5833 0 0 % confronto delle soluzioni, residui e norme [y1 y2 y3] ans = 1.0611 1.0611 1.1498 -0.0516 -0.0516 -0.1788 0.8527 0.8527 0.7480 -0.5877 -0.5877 -0.6233 0.2933 0.2933 0.5500 -0.0817 -0.0817 0 -0.4469 -0.4469 -0.6800 0.4059 0.4059 0.5833 -0.3815 -0.3815 0 0.3628 0.3628 0 norm(A'*y1-c) ans = 2.3548e-15 norm(A'*y3-c) ans = 9.0536e-16 norm(y1) ans = 1.7133 norm(y3) ans = 1.8459 % imponiamo una soluzione esatta al sistema sovradeterminato sol=ones(7,1) sol = 1 1 1 1 1 1 1 % termine noto corrispondente b=A*sol b = 2.4313 4.2575 2.1937 4.3227 3.5339 2.9243 3.3210 3.9761 2.1534 3.2807 % soluzione mediante backslash x=A\b x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 % differenza con la soluzione esatta (vettore errore) x-sol ans = 1.0e-15 * 0 0.8882 0.2220 -0.2220 -0.2220 0.2220 -0.3331 % norma dell'errore norm(x-sol) ans = 1.0474e-15 % imponiamo una soluzione esatta al sistema sottodeterminato sol=ones(10,1) sol = 1 1 1 1 1 1 1 1 1 1 % termine noto corrispondente c=A'*sol c = 5.9797 3.8909 4.9183 4.7102 5.0783 4.2131 3.6042 % soluzione mediante backslash y3=A'\c y3 = 0.2339 2.3048 1.5113 1.9847 0.4736 0 -0.2212 2.3087 0 0 % soluzione delle equazioni normali del secondo tipo z=G\c z = 0.7385 0.0580 0.8316 0.2697 -0.2919 -0.2414 0.6131 y1=A*z y1 = 0.9082 1.1807 0.9535 1.1618 0.8831 0.7021 0.8687 1.2154 0.7801 1.0591 % confronto delle soluzioni: sono diverse! [y1 y3 sol] ans = 0.9082 0.2339 1.0000 1.1807 2.3048 1.0000 0.9535 1.5113 1.0000 1.1618 1.9847 1.0000 0.8831 0.4736 1.0000 0.7021 0 1.0000 0.8687 -0.2212 1.0000 1.2154 2.3087 1.0000 0.7801 0 1.0000 1.0591 0 1.0000 % residui delle tre soluzioni norm(A'*y1-c) ans = 2.8436e-15 norm(A'*y3-c) ans = 2.2644e-15 norm(A'*sol-c) ans = 0 % norme delle tre soluzioni norm(y1) ans = 3.1165 norm(y3) ans = 4.1465 norm(sol) ans = 3.1623