CS5016: Computational Methods and Applications: Linear Systems and Interpolation

CS5016: Computational Methods and Applications
Linear Systems and Interpolation
Albert Sunny
Department of Computer Science and Engineering

Indian Institute of Technology Palakkad
01 February, 2023
Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 1 / 16

Linear systems
In applied science, one often faces equations of the following form
Ax = b
where A is a n × n square matrix whose elements are aij , x and b are
column vectors of dimension n. Component-wise, the above equation can
be written as
Xn
aij xj = bi ∀i ∈ {1, 2, . . . , n}
j=1
b can also be interpreted as a linear combination of the columns of matrix

A weighted by the vector x .

An example
Reservoir feeds water at a con-

stant pressure of 10 bar . Let
Qj = Lj ∆pj
where Lj is the length of pipe j

and is given the product of indices
of the nodes on either side of this
pipe, and ∆pj is the pressure dif-
Figure: A pipe network1 . ference across the pipe.
Write a linear system of equations
to compute the pressure at each
node.
1
“Scientific Computing with MATLAB and Octave”, Alfio Quateroni and Fausto
Saleri
Iterative solution method
An iterative method for the solution of the linear system results in a

sequence of vectors {x (k) , k ≥ 0} of Rn that converges to the exact
solution x , that is
lim x (k) = x ,
k→∞
for any given initial vector x (0) ∈ Rn .

Constructing an Iterative Method
Split matrix A, A = P − (P − A). P a suitable nonsingular matrix. Then
Px = b − (A − P )x (1)
Correspondingly, for k ≥ 0 we can define the following iterative method:
Px (k+1) = b − (A − P )x (k) (2)
(2) - (1) gives us
x (k+1) − x = (I − P −1 A)(x (k) − x ) (3)

Convergence
Let e (k) = x (k) − x denote the error at step k.

If (I − P −1 A) is a symmetric and positive definite, we have
∥e (k+1) ∥2 = ∥(I − P −1 A)e (k) ∥2 ≤ ρ(I − P −1 A)∥e (k) ∥2
where ρ(·) is known as the spectral radius (maximum modulus of

eigenvalues). If ρ(·) < 1, there is convergence.

The Jacobi method
If the diagonal entries of A are nonzero, we can set P = D , where D is

the diagonal matrix containing the diagonal entries of A. Then, we get the
following iteration
 
n
(k+1) 1 bi −
X (k)
xi = aij xj  ∀i ∈ {1, 2, . . . , n}
aii
j=1,j̸=i
Proposition
If the matrix A is strictly diagonally dominant by row, then the Jacobi
method converges.
Try proving the above proposition.

The Gauss-Siedel method
Faster convergence could be (hopefully) achieved if the new (k+1)

components already available are used
 
i−1 n
(k+1) 1  X (k+1)
X (k)
xi = bi − aij xj − aij xj  ∀i ∈ {1, 2, . . . , n}
aii
j=1 j=i+1
Caution
There are no general results stating that the Gauss-Seidel method
converges faster than Jacobi’s

Python’s numpy.linalg module
The NumPy linear algebra functions rely on BLAS and LAPACK to provide
efficient low level implementations of standard linear algebra algorithms.
To know more, visit

https://numpy.org/doc/stable/reference/routines.linalg.html

Interpolation
In several applications we may only know value of a function f at some

given points {(xi , yi ), i = 0, 1, 2, . . . , n}. How do we determine f ?
In such a situation it is natural to figure out an approximate function f˜

that satisfies the following
f˜(xi ) = yi ∀i ∈ {0, 1, 2, . . . , n}
Can you figure out a cubic function that passes through the points (0, 1),
(1, 4), (−1, 0) and (2, 15)? How many such such cubic functions exist?

A possible way!!!
Assume
f (x) = a + bx + cx 2 + dx 3
Then we get the following linear system
    
1 0 0 0 a 1
1 1 1 1   b   4 
1 −1 0 −1  c  =  0 
    
1 2 4 8 d 15
Solving the above linear system will give us coefficients of the desired
cubic polynomial.
What is the complexity of the above method?

Different kinds of interpolation
polynomial interpolation
f˜(x) = a0 + a1 x + a2 x 2 + · · · + an x n
trigonometric interpolation
f˜(x) = a−M e −iMx + · · · + a0 + · · · + aM e iMx
rational interpolation
a0 + a1 x + · · · + ak x k
f˜(x) =
b0 + b1 x + · · · + bn x n

Lagrangian polynomial interpolation
For j ∈ {0, 1, . . . , n}, define

n
Y x − xi
ψj (x) =
xj − xi
i=0,i̸=j
Note that
xj −xi
(Qn
i=0,i̸=j xj −xi = 1 if k = j
ψj (xk ) =
0 otherwise
Then, the required approximation is
n
X
f˜(x) = yj ψj (x)
j=0
Coefficients of the above polynomial can be computed in O(n2 ) time.

An example
Consider a function that passes through the points (0, 1), (1, 4), (−1, 0)
and (2, 15).
x 3 − 2x 2 − x + 2
ψ1 (x) =
2
−x 3 + x 2 + 2x
ψ2 (x) =
2
−x 3 + 3x 2 − 2x
ψ3 (x) =
6
x3 − x
ψ4 (x) =
6
and we have
f˜(x) = ψ1 (x) + 4ψ2 (x) + 15ψ4 (x)

= x3 + x2 + x + 1

Python’s scipy.interpolate module
Contains spline functions and classes, 1-D and multidimensional (univariate

and multivariate) interpolation classes, Lagrange and Taylor polynomial
interpolators, and wrappers for FITPACK and DFITPACK functions.
To know more, visit

https://docs.scipy.org/doc/scipy/reference/interpolate.html

Thank You

CS5016: Computational Methods and Applications: Linear Systems and Interpolation

Uploaded by

Copyright:

Available Formats

CS5016: Computational Methods and Applications: Linear Systems and Interpolation

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

CS5016: Computational Methods and Applications: Linear Systems and Interpolation

Uploaded by

Copyright:

Available Formats

CS5016: Computational Methods and Applications

Linear Systems and Interpolation

Department of Computer Science and Engineering

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 1 / 16

In applied science, one often faces equations of the following form

b can also be interpreted as a linear combination of the columns of matrix

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 2 / 16

Reservoir feeds water at a con-

where Lj is the length of pipe j

An iterative method for the solution of the linear system results in a

for any given initial vector x (0) ∈ Rn .

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 4 / 16

Split matrix A, A = P − (P − A). P a suitable nonsingular matrix. Then

Correspondingly, for k ≥ 0 we can define the following iterative method:

Px (k+1) = b − (A − P )x (k) (2)

(2) - (1) gives us

x (k+1) − x = (I − P −1 A)(x (k) − x ) (3)

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 5 / 16

Let e (k) = x (k) − x denote the error at step k.

∥e (k+1) ∥2 = ∥(I − P −1 A)e (k) ∥2 ≤ ρ(I − P −1 A)∥e (k) ∥2

where ρ(·) is known as the spectral radius (maximum modulus of

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 6 / 16

If the diagonal entries of A are nonzero, we can set P = D , where D is

Try proving the above proposition.

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 7 / 16

Faster convergence could be (hopefully) achieved if the new (k+1)

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 8 / 16

To know more, visit

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 9 / 16

In several applications we may only know value of a function f at some

In such a situation it is natural to figure out an approximate function f˜

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 10 / 16

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 11 / 16

f˜(x) = a−M e −iMx + · · · + a0 + · · · + aM e iMx

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 12 / 16

For j ∈ {0, 1, . . . , n}, define

Coefficients of the above polynomial can be computed in O(n2 ) time.

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 13 / 16

f˜(x) = ψ1 (x) + 4ψ2 (x) + 15ψ4 (x)

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 14 / 16

Contains spline functions and classes, 1-D and multidimensional (univariate

To know more, visit

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 15 / 16

Albert Sunny (CSE, IIT Palakkad) CS5016 01 February, 2023 16 / 16

You might also like