方法一：单位矩阵消元
1 clear; 2 n = 500; 3 A = zeros(n,n); 4 for j = 1:n 5 for i = j:n 6 A(i,j) = (i + j)^2; 7 end 8 end 9 C = A; 10 B = eye(n); 11 12 for i = 1:(n-1) 13 B(i,:) = B(i,:)/A(i,i); 14 A(i,:) = A(i,:)/A(i,i); 15 for j = (i+1):n 16 B(j,:) = B(j,:) - B(i,:)*A(j,i); 17 A(j,:) = A(j,:) - A(i,:)*A(j,i); 18 end 19 end 20 B(n,:) = B(n,:)/A(n,n); 21 A(n,:) = A(n,:)/A(n,n);
方法二：分块矩阵迭代
\left(\begin{matrix}A&C\\0&B\end{matrix}\right)^{-1}=\left(\begin{matrix}A^{-1}&-A^{-1}CB^{-1}\\0&B^{-1}\end{matrix}\right)
\left(\begin{matrix}A&0\\C&B\end{matrix}\right)^{-1}=\left(\begin{matrix}A^{-1}&0\\-B^{-1}CA^{-1}&B^{-1}\end{matrix}\right)
1 clear; 2 n = 500; 3 A = zeros(n,n); 4 for j = 1:n 5 for i = j:n 6 A(i,j) = (i + j)^2; 7 end 8 end 9 10 inv_p = 1/A(1,1); 11 for k = 2:n 12 c = A(k,1:(k-1)); 13 zero = zeros((k-1),1); 14 inv_q = [inv_p, zero; -(1/A(k,k))*c*inv_p, 1/A(k,k)]; 15 inv_p = inv_q; 16 end
第二种方法速度更快
原文：http://www.cnblogs.com/toorist/p/4818046.html