Solution of linear algebraic equations with
QR decomposition and back-substitution.
QR decomposition of a matrix A can be done via (modified)
Gram-Schmidt orthogonalization of its columns. The result is an
orthogonal matrix Q (orthogonal means QTQ=1):
do i=1,n ;
Ri,i=√(ai·ai) ; if (Ri,i=0) stop ;
qi=ai/Ri,i ← normalization of column ai
do j=i+1,n ;
Ri,j = qi·aj ; aj = aj -
qiRi,j ; enddo
← ortogonalization of all aj and qi enddo
The right triangular matrix R allows to restore the matrix A by A=QR.
The very solution of the linear system QRx=b (Q - orthogonal, R - right
triangular) is done by back-substitution
Numerical integration of two-dimensional integrals.
A two-dimensional integral of the sort
∫αβdx ∫a(x)b(x)dy f(x,y)
can be calculated by first creating a function F(x)
F(x) = ∫a(x)b(x)dy f(x,y)
using, for example, QAG, and then integrating it from α to β
using, for example, QAG. In Fortran77, though, or a copy of QAG with a
different name should be used for the second integration as recursions
are not supported.
Problems
Implement a subroutine to solve a three-diagonal linear system
QRx=b by back-substitution:
subroutine QRback(q,r,n,nmax,b)
b = QTb
do i=n,1,-1 ;
xi = ( bi - ∑(k=i+1,n) Ri,k xk ) / Ri,i ; enddo