Solve a linear system by providing `lsqr`

with a function handle that computes `A*x`

and `A'*x`

in place of the coefficient matrix `A`

.

Create an nonsymmetric tridiagonal matrix. Preview the matrix.

A = *21×21*
10 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 8 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 7 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 3 2 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 2 2 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 2 0 0 0 0 0 0 0 0 0 0
⋮

Since this tridiagonal matrix has a special structure, you can represent the operation `A*x`

with a function handle. As each row of `A`

multiplies the elements in `x`

, only a few of the results are nonzero (corresponding to the nonzeros on the tridiagonals).

The expression $\mathit{A}\text{\hspace{0.17em}}\mathit{x}$ becomes:

$\mathit{A}\text{\hspace{0.17em}}\mathit{x}=\left[\begin{array}{ccccccc}10& 2& 0& \cdots & & \cdots & 0\\ 1& 9& 2& 0& & & \vdots \\ 0& 1& \ddots & 2& 0& & \\ \vdots & 0& 1& 0& \ddots & \ddots & \vdots \\ & & 0& \ddots & 1& \ddots & 0\\ \vdots & & & \ddots & \ddots & \ddots & 2\\ 0& \cdots & & \cdots & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_{1}\\ {\mathit{x}}_{2}\\ {\mathit{x}}_{3}\\ \vdots \\ \vdots \\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}{10\mathit{x}}_{1}+2{\mathit{x}}_{2}\\ {\mathit{x}}_{1}+9{\mathit{x}}_{2}+2{\mathit{x}}_{3}\\ \vdots \\ \vdots \\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+2{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

The resulting vector can be written as the sum of three vectors:

$\mathit{A}\text{\hspace{0.17em}}\mathit{x}=\left[\begin{array}{c}{10\mathit{x}}_{1}+2{\mathit{x}}_{2}\\ {\mathit{x}}_{1}+9{\mathit{x}}_{2}+2{\mathit{x}}_{3}\\ \vdots \\ \vdots \\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+2{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_{1}\\ {\mathit{x}}_{2}\\ \vdots \\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}{10\mathit{x}}_{1}\\ {9\mathit{x}}_{2}\\ \vdots \\ 9{\mathit{x}}_{20}\\ 10{\mathit{x}}_{21}\end{array}\right]+2\cdot \left[\begin{array}{c}{\mathit{x}}_{2}\\ {\mathit{x}}_{3}\\ \vdots \\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

Likewise, the expression for ${\mathit{A}}^{\mathit{T}}\text{\hspace{0.17em}}\mathit{x}$ becomes:

${\mathit{A}}^{\mathit{T}}\text{\hspace{0.17em}}\mathit{x}=\left[\begin{array}{ccccccc}10& 1& 0& \cdots & & \cdots & 0\\ 2& 9& 1& 0& & & \vdots \\ 0& 2& \ddots & 1& 0& & \\ \vdots & 0& 2& 0& \ddots & \ddots & \vdots \\ & & 0& \ddots & 1& \ddots & 0\\ \vdots & & & \ddots & \ddots & \ddots & 1\\ 0& \cdots & & \cdots & 0& 2& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_{1}\\ {\mathit{x}}_{2}\\ {\mathit{x}}_{3}\\ \vdots \\ \vdots \\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}{10\mathit{x}}_{1}+{\mathit{x}}_{2}\\ {2\mathit{x}}_{1}+9{\mathit{x}}_{2}+{\mathit{x}}_{3}\\ \vdots \\ \vdots \\ {2\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {2\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

${\mathit{A}}^{\mathit{T}}\text{\hspace{0.17em}}\mathit{x}=\left[\begin{array}{c}{10\mathit{x}}_{1}+{\mathit{x}}_{2}\\ {2\mathit{x}}_{1}+9{\mathit{x}}_{2}+{\mathit{x}}_{3}\\ \vdots \\ \vdots \\ {2\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {2\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]=2\cdot \left[\begin{array}{c}0\\ {\mathit{x}}_{1}\\ {\mathit{x}}_{2}\\ \vdots \\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}{10\mathit{x}}_{1}\\ {9\mathit{x}}_{2}\\ \vdots \\ 9{\mathit{x}}_{20}\\ 10{\mathit{x}}_{21}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_{2}\\ {\mathit{x}}_{3}\\ \vdots \\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

In MATLAB®, write a function that creates these vectors and adds them together, thus giving the value of `A*x`

or `A'*x`

, depending on the flag input:

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
y = [0; x(1:20)] ...
+ [(10:-1:0)'; (1:10)'].*x ...
+ 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
y = 2*[0; x(1:20)] ...
+ [(10:-1:0)'; (1:10)'].*x ...
+ [x(2:end); 0];
end
end

(This function is saved as a local function at the end of the example.)

Now, solve the linear system $\mathrm{Ax}=\mathit{b}$ by providing `lsqr`

with the function handle that calculates `A*x`

and `A'*x`

. Use a tolerance of `1e-6`

and 25 iterations. Specify $\mathit{b}$ as the row sums of $\mathit{A}$ so that the true solution for $\mathit{x}$ is a vector of ones.

lsqr converged at iteration 21 to a solution with relative residual 5.4e-13.

x1 = *21×1*
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
⋮