Linear Algebra
|
This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production. The reference material consulted while preparing these pages is listed in the section bibliography. |
Written against MATLAB R2026b. The snippets in the "Octave 11.3.0 output" blocks were run in GNU Octave 11.3.0 using only code that Octave also accepts; the MATLAB layout of the same numbers may differ in spacing and in the last digits of round-off values. The decomposition object and the exact MATLAB warning texts are not available in Octave, so those parts are checked against the R2026b reference pages and labelled "not run here". Everything on this page is base MATLAB: no toolbox is needed. Functions of the Statistics, Optimization and Symbolic Math toolboxes that overlap with this topic are listed on Toolboxes and Add-On Products.
MATLAB began as an interactive front end to the LINPACK and EISPACK linear algebra libraries (see the introduction), and linear algebra is still the centre of the language. This page is about the dense, full-storage case. Sparse storage and iterative solvers are on Sparse Matrices and Graphs; the array basics (creation, indexing, operators) are on Matrices and Arrays and Operators and Precedence; the floating-point background of every accuracy statement below is on Numeric Types and Floating Point.
Solving A*x = b: backslash and slash
To solve the system write x = A\b. The operator is called left division (mldivide). Its mirror image x = b/A (mrdivide) solves , and the two are related by b/A = (A'\b')'. The documentation explains the idea with the scalar equation 7x = 21: you divide, you do not compute the reciprocal of 7 and multiply. MATLAB solves matrix systems the same way, without forming the inverse.
The coefficient matrix may have any shape, and the same operator covers three cases:
Shape of A (m-by-n) |
What A\b returns |
|---|---|
Square (m = n), nonsingular |
The unique solution of the system. |
Overdetermined (m > n) |
The least-squares solution, the |
Underdetermined (m < n) |
A basic solution with at most m nonzero components. For the minimum-norm solution use |
A small, well-conditioned example (a symmetric tridiagonal matrix, the shape that appears in the plate section below). The right-hand side was chosen so that the exact answer is a vector of ones:
A = [4 -1 0; -1 4 -1; 0 -1 4];
b = [3; 2; 3];
x = A\b
r = b - A*x
c = [1 2 3];
y = c/A
norm(y' - A'\c')Octave 11.3.0 output:
x =
1
1
1
r =
0
8.8818e-16
0
y =
0.4643 0.8571 0.9643
ans = 0The residual r is not exactly zero: it is round-off at the level of eps. The last line confirms the identity between c/A and the transposed backslash.
Why not inv(A)*b
The inverse is almost never needed. The documentation says that computing it is more work and, in finite-digit arithmetic, less accurate, and the inv reference page recommends x = A\b over x = inv(A)*b. The reference page for mldivide adds that for a square matrix A\B is roughly equal to inv(A)*B but is processed differently and more robustly.
The Hilbert matrix is the classic ill-conditioned test case. Building b = A*ones(10,1) makes the exact answer known, so the error of each method can be measured:
n = 10;
A = hilb(n);
xtrue = ones(n,1);
b = A*xtrue;
fprintf("cond(A) = %.3g\n", cond(A));
fprintf("backslash error = %.3g\n", norm(A\b - xtrue));
fprintf("inv(A)*b error = %.3g\n", norm(inv(A)*b - xtrue));
fprintf("backslash residual = %.3g\n", norm(b - A*(A\b)));
fprintf("inv(A)*b residual = %.3g\n", norm(b - A*(inv(A)*b)));Octave 11.3.0 output:
cond(A) = 1.6e+13
backslash error = 0.00212
inv(A)*b error = 0.00333
backslash residual = 6.08e-16
inv(A)*b residual = 0.000371Both errors are large, because the condition number is about and no algorithm can recover digits the data never carried. The difference shows in the residual: the backslash solution reproduces b to round-off, while the inverse-based one leaves a residual nine orders of magnitude larger. The exact numbers depend on the BLAS library of the platform and will differ a little in MATLAB. The documentation also warns that a small residual does not guarantee an accurate solution when the matrix is ill conditioned, which is the case here for the backslash error.
Two further reasons not to form the inverse: it costs more operations than a factorization, and for structured matrices (triangular, symmetric positive definite, banded, sparse) the inverse destroys the structure that mldivide exploits. When you genuinely need the inverse (for example the covariance of an estimator), inv exists and warns when the matrix is badly scaled or nearly singular.
Source: MATLAB documentation: inv, mldivide
How mldivide chooses a solver
mldivide is not one algorithm. The documentation says it chooses a solver based on the properties of A (and of B) so as to exploit symmetry or other structure, and that the choice is shown as a flowchart on the reference page. The page does not describe the branches in text, so the flowchart below is an original simplified illustration of the ordering of the checks, not a copy of the MathWorks figure; the solver families agree with the types accepted by the decomposition object (diagonal, triangular, permutedTriangular, banded, hessenberg, chol, ldl, lu, qr, cod). Read the reference page for the exact conditions.
The practical consequences:
-
Structure pays. The same
A\bis a few multiplications for a triangular matrix, about one third of the work of LU for a symmetric positive definite one, and much cheaper again for a banded or sparse one. Keep structure visible: build tridiagonal matrices withspdiagsrather than a full array when the size is large. -
Square systems that are singular or nearly so produce a warning, because the solver estimates the reciprocal condition number
RCONDas a by-product. -
Rank-deficient non-square systems are solved in the least-squares sense and return a basic solution (see the
null,orthandpinvsection).
Decomposition objects
When you solve many systems with the same coefficient matrix, factor once and reuse the factors. The decomposition function returns an object that stores a factorization and solves with it:
dA = decomposition(A); % type chosen automatically
x1 = dA\b1;
x2 = dA\b2; % no refactorization
cond_ill = isIllConditioned(dA); % uses the stored information
dS = decomposition(A,"chol"); % force a typeNot run here (Octave 11.3.0 has no decomposition). The reference page lists the types "auto" (default), "qr", "cod", "lu", "ldl", "chol", "triangular", "permutedTriangular", "banded", "hessenberg" and "diagonal". Solving is done with the ordinary operators: dA\b and b/dA, so existing code needs only the first line changed. The object has the read-only properties MatrixSize, Type, Datatype, IsConjugateTransposed, IsReal, IsSparse and ScaleFactor, a settable CheckCondition, and the methods isIllConditioned, rcond and (for the "qr" and "cod" types) rank. The page states that it is typically much faster than A\b for repeated solves and shows an inverse-iteration example whose run time drops from about 8.9 s to about 1.0 s. It also lets you pick the type by hand, for instance "qr" to obtain a least-squares solution where the automatic choice would give NaN values. The object is most valuable inside loops, such as the time-stepping of an implicit heat equation scheme (see Integration, ODEs and Optimization).
Determinant, rank, trace and the singular-matrix warnings
A square matrix is singular when it has no inverse. Four functions look at that question from different angles, and they do not agree about what to trust:
| Function | What it tells you |
|---|---|
|
The determinant, computed from the LU factors. The reference page says to avoid using it to check singularity: the determinant can be arbitrarily close to zero without conveying anything about singularity. Scaling a matrix by a constant scales the determinant by that constant to the n-th power. |
|
The number of singular values above a tolerance, computed with the SVD. The default tolerance is |
|
An inexpensive estimate of the reciprocal of the 1-norm condition number. Near 1 means well conditioned, near 0 means close to singular. The reference page recommends |
|
The condition number for inversion, the ratio of the largest to the smallest singular value in the default 2-norm. It is 1 for a perfectly conditioned matrix and |
The classic demonstration pairs magic(4) (singular in exact arithmetic, but its floating-point det is a tiny nonzero number) with a scaled identity (its det is tiny, yet it is perfectly conditioned):
A = magic(4);
fprintf("det = %.4g\n", det(A));
fprintf("rank = %d\n", rank(A));
fprintf("rcond = %.3g\n", rcond(A));
d = det(0.1*eye(100));
fprintf("det(0.1*eye(100)) = %.3g, rank = %d, rcond = %g\n", d, rank(0.1*eye(100)), rcond(0.1*eye(100)));
B = [1 0; 0 0]; b = [1; 1];
x = B\bOctave 11.3.0 output (the warning text is Octave’s, with its "called from" lines removed; MATLAB prints its own wording, below):
det = -1.45e-12
rank = 3
rcond = 1.31e-17
det(0.1*eye(100)) = 1e-100, rank = 100, rcond = 1
warning: matrix singular to machine precision
x =
1
0det(magic(4)) is about , not zero, while the matrix has rank 3 and an rcond of order . In the other direction det(0.1*eye(100)) is , which looks singular, yet rcond is exactly 1. The determinant answers the wrong question. The third matrix is exactly singular: B\b returns a result of some kind, but the system has no exact solution, so the warning matters more than the printed values.
The MATLAB warnings. The mldivide reference page describes two cases, both reported with a warning and the value of RCOND:
| Condition | Message (from the reference page) | Warning identifier (not run here) |
|---|---|---|
|
"Matrix is close to singular or badly scaled. Results may be inaccurate." followed by |
|
|
"Matrix is singular to working precision."; the result can contain |
|
[msg,id] = lastwarn right after the warning returns the message and the identifier, which is how to confirm which of the two was raised. Because they are warnings and not errors, a script that must not continue with an unreliable answer can turn them into one:
warning("error","MATLAB:singularMatrix"); % not run here: makes the warning throw
try
x = B\b;
catch err
disp(err.identifier)
endThe trace function sums the diagonal, and norm measures size: for a vector norm(v) is the Euclidean length, norm(v,1) the sum of absolute values and norm(v,Inf) the largest absolute value; for a matrix norm(A) is the largest singular value, norm(A,1) the maximum absolute column sum, norm(A,Inf) the maximum absolute row sum and norm(A,"fro") the Frobenius norm. For sparse matrices use norm(A,"fro"), because the 2-norm is not supported for sparse input.
A = [1 2; 3 4];
v = [3 4];
fprintf("trace(A) = %d\n", trace(A));
fprintf("norm(v) = %g, norm(v,1) = %g, norm(v,Inf) = %g\n", norm(v), norm(v,1), norm(v,Inf));
fprintf("norm(A) = %.4f, norm(A,1) = %g, norm(A,'fro') = %.4f\n", norm(A), norm(A,1), norm(A,'fro'));
fprintf("cond(A) = %.4f, rcond(A) = %.4f\n", cond(A), rcond(A));
H = hilb(8);
fprintf("cond(hilb(8)) = %.3g, rcond(hilb(8)) = %.3g\n", cond(H), rcond(H));Octave 11.3.0 output:
trace(A) = 5
norm(v) = 5, norm(v,1) = 7, norm(v,Inf) = 4
norm(A) = 5.4650, norm(A,1) = 6, norm(A,'fro') = 5.4772
cond(A) = 14.9330, rcond(A) = 0.0476
cond(hilb(8)) = 1.53e+10, rcond(hilb(8)) = 2.95e-11cond and rcond move in opposite directions and are not equal; rcond is cheaper but, in the words of the reference page, less reliable. A rule of thumb from the numerical analysis the documentation relies on: expect to lose about decimal digits of accuracy when solving with A.
Null space, range and the pseudoinverse
When A does not have full rank, three SVD-based functions describe what is missing and what remains:
-
Z = null(A)returns an orthonormal basis of the null space: the columns ofVthat belong to zero singular values.null(A,"rational")returns a basis from the reduced row echelon form, usually not orthonormal and less accurate (the older'r'spelling still works). -
Q = orth(A)returns an orthonormal basis of the range; it hasrank(A)columns. -
pinv(A)is the Moore-Penrose pseudoinverse, defined for every matrix and computed fromsvd(A). Singular values at or belowmax(size(A))*eps(norm(A))are treated as zero (the same default tolerance asrank,nullandorth). For a square nonsingular matrixpinv(A)is "simply an expensive way" to computeinv(A), in the words of its reference page.
All the general solutions of a consistent system are one particular solution plus anything in the null space. For an inconsistent system pinv(A)*b gives the least-squares solution of minimum norm. The page also points to lsqminnorm(A,b) as the more efficient way to get the same minimum-norm solution.
A = [1 2 3; 4 5 6; 7 8 9];
fprintf("rank = %d\n", rank(A));
Z = null(A);
fprintf("size(null) = %s, norm(A*Z) = %.2e\n", mat2str(size(Z)), norm(A*Z));
Q = orth(A);
fprintf("size(orth) = %s, orthonormal: %.2e\n", mat2str(size(Q)), norm(Q'*Q - eye(2)));
b = [1; 2; 3];
x = pinv(A)*b;
fprintf("pinv solution = %s, residual = %.2e\n", mat2str(x',4), norm(A*x - b));
R = [6 8 7 3; 3 5 4 1]; c = [7; 8];
p = R\c;
fprintf("R\\c = %s, norm = %.4f, nnz = %d\n", mat2str(p',4), norm(p), nnz(abs(p)>1e-12));
q = pinv(R)*c;
fprintf("pinv = %s, norm = %.4f\n", mat2str(q',4), norm(q));
fprintf("both solve: %.2e %.2e\n", norm(R*p-c), norm(R*q-c));Octave 11.3.0 output:
rank = 2
size(null) = [3 1], norm(A*Z) = 4.44e-16
size(orth) = [3 2], orthonormal: 5.33e-16
pinv solution = [-0.05556 0.1111 0.2778], residual = 2.40e-15
R\c = [-1.511 2.664 0.5766 -3.095], norm = 4.3923, nnz = 4
pinv = [-1.511 2.664 0.5766 -3.095], norm = 4.3923
both solve: 7.94e-15 7.16e-15The matrix [1 2 3; 4 5 6; 7 8 9] has rank 2 (a null space of dimension one and a range of dimension two), and b = [1;2;3] lies in its range, so pinv(A)*b solves the system to round-off. For the underdetermined 2-by-4 system R this particular Octave build returned the same vector from R\c and from pinv(R)*c. That is not guaranteed in general: the documentation describes R\c as a basic solution with at most rank® nonzero components and pinv(R)*c as the minimum-norm solution, and they are different vectors unless the basic solution happens to be the minimum-norm one. Do not rely on them being equal. In this run nnz counts four nonzero entries in p, so it was not a basic solution at all, which shows that the behaviour of the underdetermined case is implementation specific. Use lsqminnorm when you need the minimum-norm answer and a documented guarantee.
Source: MATLAB documentation: null, orth, pinv, lsqminnorm
The factorizations: lu, qr and chol
Three factorizations underlie most of linear algebra, and mldivide uses them internally. You call them directly to reuse factors, to inspect pivoting or to implement your own algorithms.
| Function | Factorization | Solving |
|---|---|---|
|
with |
|
|
with |
Least squares: |
|
with |
|
The factoring cost is proportional to , and each additional solve with the stored factors costs . That is what a decomposition object automates.
A = [2 1 1; 4 3 3; 8 7 9];
[L,U,P] = lu(A)
fprintf("norm(P*A - L*U) = %.2e\n", norm(P*A - L*U));
[Q,R] = qr(A);
fprintf("norm(A - Q*R) = %.2e, norm(Q'*Q - I) = %.2e\n", norm(A-Q*R), norm(Q'*Q-eye(3)));
S = [4 -1 0; -1 4 -1; 0 -1 4];
R = chol(S)
fprintf("norm(R'*R - S) = %.2e\n", norm(R'*R - S));
N = [1 2; 2 1];
[R2,flag] = chol(N);
fprintf("chol(N) second output flag = %d\n", flag);
E = [1 1; 1 2; 1 3; 1 4]; y = [1;2;2;4];
[Q1,R1] = qr(E,0);
c = R1\(Q1'*y);
fprintf("least squares via qr: %s ; via backslash: %s\n", mat2str(c',6), mat2str((E\y)',6));Octave 11.3.0 output:
L =
1.0000 0 0
0.2500 1.0000 0
0.5000 0.6667 1.0000
U =
8.0000 7.0000 9.0000
0 -0.7500 -1.2500
0 0 -0.6667
P =
Permutation Matrix
0 0 1
1 0 0
0 1 0
norm(P*A - L*U) = 0.00e+00
norm(A - Q*R) = 4.83e-15, norm(Q'*Q - I) = 3.04e-16
R =
2.0000 -0.5000 0
0 1.9365 -0.5164
0 0 1.9322
norm(R'*R - S) = 1.11e-16
chol(N) second output flag = 2
least squares via qr: [-4.44089e-16 0.9] ; via backslash: [7.54952e-16 0.9]The pivot of the first column is 8, the largest, so P moves the third row to the top. The matrix N = [1 2; 2 1] is symmetric but indefinite (its eigenvalues are 3 and -1), so chol stops at the second pivot and reports flag = 2. chol can therefore be used as a cheap positive-definiteness test. The straight-line fit through E and y gives the same intercept (zero to round-off) and slope 0.9 from the QR route and from backslash, which uses a QR-based solver for non-square systems. In MATLAB, P from lu is a plain double matrix (Octave prints it as "Permutation Matrix" and the layout of everything above can differ).
Source: MATLAB documentation: Factorizations, lu, qr, chol
Eigenvalues and singular values
For a square matrix, lambda = eig(A) returns the eigenvalues as a column vector and [V,D] = eig(A) returns the eigenvectors as the columns of V and the eigenvalues on the diagonal of D, so that A*V = V*D. The eigenvectors are scaled to unit length. Some matrices (defective ones) do not have a full set of eigenvectors; for those, the documentation recommends the Schur decomposition [U,S] = schur(A), which is more robust.
The generalised problem is eig(A,B), and the right eigenvectors satisfy A*V = B*V*D. Further options on the reference page: "balance" (default) and "nobalance" for the preliminary balancing step (the page recommends "nobalance" when the entries differ greatly in scale), "vector" or "matrix" for the form of the eigenvalues, "chol" or "qz" for the algorithm of the generalised problem, and a third output W with the left eigenvectors.
The singular value decomposition [U,S,V] = svd(A) factors any m-by-n matrix as with unitary U and V and a diagonal S of non-negative singular values in decreasing order; s = svd(A) returns just the singular values. When A is tall, svd(A,"econ") drops the columns of U that would only multiply zeros: U is m-by-n and S is n-by-n. (The form svd(A,0) is the older spelling and is not recommended.) Singular values are the bridge to the earlier section: rank counts those above a tolerance, cond is the ratio of the largest to the smallest, norm(A) is the largest, and pinv, null and orth are built from them. For large sparse matrices svds computes a few of them, and svdsketch builds a low-rank approximation.
A = [2 0; 1 3];
[V,D] = eig(A)
S = [2 -1 0; -1 2 -1; 0 -1 2];
e = eig(S);
fprintf("%.4f ", e); fprintf("\n");
B = diag([1 2 3]);
[V2,D2] = eig(S,B);
fprintf("generalized eig(S,B): %s\n", mat2str(diag(D2)',4));
fprintf("check A*V = B*V*D: %.2e\n", norm(S*V2 - B*V2*D2));
M = [1 2; 3 4; 5 6];
s = svd(M)
[U,Sg,W] = svd(M);
fprintf("size U %s S %s V %s, norm(M-U*S*V') = %.2e\n", mat2str(size(U)), mat2str(size(Sg)), mat2str(size(W)), norm(M-U*Sg*W'));
[U1,S1,W1] = svd(M,"econ");
fprintf("econ sizes U %s S %s V %s\n", mat2str(size(U1)), mat2str(size(S1)), mat2str(size(W1)));
fprintf("cond(M) = %.4f = s(1)/s(end) = %.4f\n", cond(M), s(1)/s(end));
fprintf("rank by counting s > tol: %d\n", sum(s > max(size(M))*eps(s(1))));Octave 11.3.0 output:
V =
0 0.7071
1.0000 -0.7071
D =
Diagonal Matrix
3 0
0 2
0.5858 2.0000 3.4142
generalized eig(S,B): [0.2792 1 2.387]
check A*V = B*V*D: 5.68e-16
s =
9.5255
0.5143
size U [3 3] S [3 2] V [2 2], norm(M-U*S*V') = 1.30e-15
econ sizes U [3 2] S [2 2] V [2 2]
cond(M) = 18.5213 = s(1)/s(end) = 18.5213
rank by counting s > tol: 2The first matrix is triangular, so its eigenvalues are its diagonal entries 2 and 3, but the order in which eig lists them is not specified, so do not rely on it (Octave lists 3 before 2 here for D). The matrix S is the 1-D second-difference matrix, with eigenvalues : 0.5858, 2 and 3.4142. The eigenvalues of this kind of matrix reappear in the plate below.
Source: MATLAB documentation: Eigenvalues, Singular Values, eig, svd
Matrix functions: expm, logm, sqrtm and kron
Functions that end in m act on the matrix as a whole; their element-wise cousins act on each entry. The two are not interchangeable.
-
expm(A)is the matrix exponential, , and it solves the linear differential system as . It is computed by scaling and squaring with a Padé approximation, not via eigenvectors.exp(A)is the element-wise exponential and gives a different matrix. -
logm(A)is the principal matrix logarithm, the inverse ofexpm. IfAis singular or has eigenvalues on the negative real axis the principal logarithm is undefined, andlogmreturns a non-principal one and warns. Some matrices have no logarithm at all, such as[0 1; 0 0]. -
sqrtm(A)is the principal square root, the matrixXwithX*X = A. The result is complex ifAhas eigenvalues with negative real parts. Element-wisesqrt(A)is a different matrix. -
kron(A,B)is the Kronecker product, the block matrix whose block (i,j) isA(i,j)*B; its size ism*p-by-n*qfor an m-by-n and a p-by-q input. If either input is sparse, the result is sparse. It is the tool for building discretised operators in several dimensions, as the plate shows next.
A = [0 1; -1 0];
E = expm(A)
fprintf("expm(A) vs rotation: %.2e\n", norm(E - [cos(1) sin(1); -sin(1) cos(1)]));
fprintf("exp(A) elementwise: %s\n", mat2str(exp(A),4));
L = logm(E);
fprintf("logm(expm(A)) - A: %.2e\n", norm(L - A));
S = [4 1; 1 3];
R = sqrtm(S);
fprintf("norm(R*R - S) = %.2e\n", norm(R*R - S));
fprintf("elementwise sqrt differs: %s\n", mat2str(sqrt(S),4));
K = kron([1 2; 3 4], eye(2))
J = kron(eye(2), [1 2; 3 4]);
disp(J)Octave 11.3.0 output:
E =
0.5403 0.8415
-0.8415 0.5403
expm(A) vs rotation: 1.11e-16
exp(A) elementwise: [1 2.718;0.3679 1]
logm(expm(A)) - A: 3.50e-16
norm(R*R - S) = 2.44e-15
elementwise sqrt differs: [2 1;1 1.732]
K =
1 0 2 0
0 1 0 2
3 0 4 0
0 3 0 4
1 2 0 0
3 4 0 0
0 0 1 2
0 0 3 4The generator [0 1; -1 0] is a rotation, so expm(A) is the rotation matrix by one radian, while the element-wise exp(A) has no such meaning. kron(X,eye(2)) interleaves the entries of X; kron(eye(2),X) puts copies of X on the diagonal. These are the two building blocks of the next section.
The plate’s steady state as A\b
The heatflow plate reaches a steady state when its temperature stops changing. For heat diffusion this means the Laplace equation in the interior, with the temperature fixed on the edges (the Dirichlet boundary kind of heatflow.BoundaryKind). On a uniform grid with equal spacing in both directions, the five-point approximation turns it into one equation per interior point:
where a neighbour that lies on the boundary is a known number and moves to the right-hand side. For an n-by-n interior grid, unknowns taken column by column, the coefficient matrix is
which is exactly kron(I,T1) + kron(T1,I). The matrix is n^2-by-n^2, symmetric, positive definite and has at most five entries per row, so mldivide sends it down the Cholesky path (or the sparse Cholesky path when it is stored as sparse, as spdiags and speye make it). The right-hand side is zero everywhere except at interior points adjacent to a hot edge, where it equals the boundary temperature. Take the top edge at 100 degrees Celsius and the other three at 0:
% steady state of a plate: interior n-by-n grid, top edge 100, other edges 0
n = 4;
e = ones(n,1);
T1 = spdiags([-e 2*e -e], -1:1, n, n);
I = speye(n);
A = kron(I,T1) + kron(T1,I);
fprintf("size(A) = %s, nnz = %d, symmetric = %d\n", mat2str(size(A)), nnz(A), issymmetric(A));
top = 100;
b = zeros(n,n);
b(1,:) = top;
b = b(:);
t = A\b;
Tin = reshape(t,n,n)
fprintf("residual = %.2e\n", norm(A*t-b));
fprintf("chol succeeds: %d\n", ~isempty(chol(full(A))));
fprintf("cond(full(A)) = %.3f\n", cond(full(A)));
fprintf("eig min/max: %.4f %.4f\n", min(eig(full(A))), max(eig(full(A))));
n = 40;
e = ones(n,1); T1 = spdiags([-e 2*e -e], -1:1, n, n); I = speye(n);
A = kron(I,T1) + kron(T1,I);
b = zeros(n,n); b(1,:) = 100; b = b(:);
t = A\b; Tin = reshape(t,n,n);
fprintf("n=40: unknowns %d, T(n/2,n/2) = %.4f, max interior %.4f\n", n*n, Tin(n/2,n/2), max(t));Octave 11.3.0 output:
size(A) = [16 16], nnz = 64, symmetric = 1
Tin =
45.4545 59.4697 59.4697 45.4545
22.3485 32.9545 32.9545 22.3485
10.9848 17.0455 17.0455 10.9848
4.5455 7.1970 7.1970 4.5455
residual = 1.12e-13
chol succeeds: 1
cond(full(A)) = 9.472
eig min/max: 0.7639 7.2361
n=40: unknowns 1600, T(n/2,n/2) = 26.0169, max interior 95.0818Row 1 of Tin is the row next to the hot edge, and the temperature falls monotonically away from it and is symmetric left to right, as the physics requires. The matrix is well behaved: it is positive definite (chol succeeds, and the smallest eigenvalue is 0.7639), and its condition number is 9.47 for the 16 unknowns. That number grows like the square of the number of grid points per side as the grid is refined, so the 1600-unknown system is still solved accurately by A\b but is the place to switch to the iterative solvers of Sparse Matrices and Graphs when the grid becomes very fine. In a project, this computation would live in a function of the heatflow namespace that builds the sparse matrix with kron and returns the temperature field for a heatflow.Plate.
Source: MATLAB documentation: kron, mldivide
Where older books differ
The books this section follows were written for much older releases, and on this topic the differences are mostly in what they omit. Neither book describes decomposition objects (R2017b), which replace the habit of calling lu by hand and keeping the factors in variables for repeated solves. Older texts show qr(A,0) and svd(A,0) for the economy size; both still work, but the reference pages now call them not recommended and prefer the "econ" option. Older texts show null(A,'r'), which was renamed "rational" in R2022a (the old spelling still works). Examples that compute inv(A)*b or use det(A) == 0 to test singularity are common in older material and should be replaced with A\b and with rcond or cond. The pinv page now recommends lsqminnorm for solving. Finally, the behaviour of underdetermined A\b (a basic solution, not the minimum-norm one) is a frequent source of confusion in older books; the current documentation states it explicitly.