ASSIGNMENT 2¶

In [1]:
import pandas as pd
import numpy as np
import math

Question 1¶

Solve the equations [ \begin{cases} 2x - y + z = 4\\[4pt] x + y + z = 3\\[4pt] 3x - y - z = 1 \end{cases} ]

using the MATLAB’s left-division operator (/). Check your solution by computing the residual. Also compute the determinant (det) and the condition estimator (rcond). What do you conclude?

In [8]:
# A: the coefficient matrix
# b: the constant vector
A = np.array([[2, -1,  1],
              [1,  1,  1],
              [3, -1, -1]], dtype=float)
b = np.array([4, 3, 1], dtype=float)

# Ax = b
x = np.linalg.solve(A, b)

# compute residual r and its norm
r = A @ x - b
residual_norm = np.linalg.norm(r)

# compute determinant and condition estimator
det_A = np.linalg.det(A)
cond_A = np.linalg.cond(A, 2)
rcond_A = 1 / cond_A  # reciprocal condition estimate

print("Solution vector x:\n", x)
print("\nResidual vector r:\n", r)
print("\nResidual norm ||r|| =\n", residual_norm)
print("\nDeterminant det(A) =", det_A)
print("Condition number cond(A) =", cond_A)
print("Reciprocal condition estimate rcond(A) =", rcond_A)

Solution vector x:
 [1. 0. 2.]

Residual vector r:
 [0. 0. 0.]

Residual norm ||r|| =
 0.0

Determinant det(A) = -7.999999999999998
Condition number cond(A) = 3.577608106737146
Reciprocal condition estimate rcond(A) = 0.27951636125736007
In [9]:
# comments based on the results
if rcond_A > 1e-3:
    print("\nThe system is well-conditioned; the solution is reliable.")
else:
    print("\nThe system is ill-conditioned; numerical errors may be significant.")

The system is well-conditioned; the solution is reliable.

Question 2¶

This problem, suggested by R.V. Andree, demonstrates ill-conditioning (where small changes in the coefficients cause large changes in the solution). Use the MATLAB left division operator ( / ) to show that the solution of the system:

\begin{cases} x + 5.000\,y = 17.0\\[2pt] 1.5\,x + 7.501\,y = 25.503 \end{cases}

Compute the residual.

Then change the right-hand side of to 25.501, 25.502, 25.504, etc. and compute the new solutions and residuals.

If the coefficients are subject to experimental errors, the solution is clearly meaningless. Use rcond to find the condition estimator and det to compute the determinant. Do these values confirm ill conditioning?

Another way to anticipate ill conditioning is to perform a sensitivity analysis on the coefficients: change them all in turn by the same small percentage, and observe what effect this has on the solution.

In [18]:
A = np.array([[1.0, 5.000],
              [1.5, 7.501]], dtype=float)
b = np.array([17.0, 25.503], dtype=float)

x = np.linalg.solve(A, b)
r = A @ x - b

print("Original solution (x, y):", x)
print("Residual:", r)

Original solution (x, y): [2. 3.]
Residual: [0. 0.]
In [19]:
b_changes = [25.501, 25.502, 25.504]
for b_changed in b_changes:
    b_new = np.array([17.0, b_changed])
    x_new = np.linalg.solve(A, b_new)
    r_new = A @ x_new - b_new
    print(f"\nRight-hand side = {b_new}")
    print("Solution (x, y):", x_new)
    print("Residual:", r_new)

Right-hand side = [17.    25.501]
Solution (x, y): [12.  1.]
Residual: [ 0.00000000e+00 -3.55271368e-15]

Right-hand side = [17.    25.502]
Solution (x, y): [7. 2.]
Residual: [0. 0.]

Right-hand side = [17.    25.504]
Solution (x, y): [-3.  4.]
Residual: [0. 0.]
In [20]:
det_A = np.linalg.det(A)
cond_A = np.linalg.cond(A, 2)
rcond_A = 1 / cond_A

print("\nDeterminant det(A) =", det_A)
print("Condition number cond(A) =", cond_A)
print("Reciprocal condition estimate rcond(A) =", rcond_A)

Determinant det(A) = 0.00099999999999989
Condition number cond(A) = 84515.00098783901
Reciprocal condition estimate rcond(A) = 1.1832218994399485e-05
In [38]:
# sensitivity analysis: change each coefficient by +0.1%
change = 0.001
print("\nSensitivity Analysis")
print("\nSolutions after increasing the values of A by 0.1%:\n")
for i in range(A.shape[0]):
    for j in range(A.shape[1]):
        A_changed = A.copy()
        A_changed[i,j] *= (1 + change)
        x_changed = np.linalg.solve(A_changed, b)
        print(f"The new value of A[{i},{j}]={round(A_changed[i,j],3)} -> solution: {x_changed}")

Sensitivity Analysis

Solutions after increasing the values of A by 0.1%:

The new value of A[0,0]=1.001 -> solution: [0.23526644 3.35289966]
The new value of A[0,1]=5.005 -> solution: [19.31       -0.46153846]
The new value of A[1,0]=1.501 -> solution: [-0.30769231  3.46153846]
The new value of A[1,1]=7.509 -> solution: [15.23550171  0.35289966]

Results:¶

  1. Original system: The solution is $(x,y) = (2,3)$ with a residual close to zero, which confirms that the solution satisfies the equations.

  2. Small changes in the right-hand side:

    • Changing the second equation's right-hand side from 25.503 → 25.501 produces a noticeably different solution.
    • Further small changes (25.502, 25.504) also change the solution significantly.
    • This shows the system is extremely sensitive to tiny perturbations, and the system is ill-conditioned.

  3. Determinant and condition number:

    • det(A) is small, indicating the matrix is near singular.
    • cond(A) is very large and rcond is very small (near 0), confirming that the system is ill-conditioned.

  4. Sensitivity analysis on coefficients:

    • Increasing any single coefficient by $0.1%$ leads to large variations in the solution, further demonstrating the system’s sensitivity.

Conclusion: The system is ill-conditioned. Even very small changes or errors in coefficients or measurements can produce completely unreliable solutions.

Question 3¶

If you are familiar with Gauss reduction it is an excellent programming exercise to code a Gauss reduction directly with operations on the rows of the augmented coefficient matrix. See if you can write a function

x = mygauss(a, b)

to solve the general system Ax = b. Skillful use of the colon operator in the row operations can reduce the code to a few lines! Test it on A and b with random entries, and on the systems in Questions 1 and 2.

In [40]:
def mygauss(A, b):
    """
    solve Ax = b
    A: the coefficient matrix of size n x n
    b: the right-hand-side vector of size n
    x: the solution vector of size n
    """
    A = A.astype(float).copy()
    b = b.astype(float).copy()
    n = len(b)

    # forward elimination
    for i in range(n):
        # partial pivoting: swap row with the largest pivot
        max_row = np.argmax(np.abs(A[i:, i])) + i
        if max_row != i:
            A[[i, max_row]] = A[[max_row, i]]
            b[[i, max_row]] = b[[max_row, i]]

        # eliminate entries below pivot
        for j in range(i + 1, n):
            factor = A[j, i] / A[i, i]
            A[j, i:] -= factor * A[i, i:]
            b[j] -= factor * b[i]

    # back substitution
    x = np.zeros(n)
    for i in reversed(range(n)):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]
    return x
In [50]:
# test on random matrix and vectors
np.random.seed(522)
A_rand = np.random.rand(3,3)
b_rand = np.random.rand(3)
x_rand = mygauss(A_rand, b_rand)
print("Solution x =", x_rand)
print("Residual =", np.linalg.norm(A_rand @ x_rand - b_rand))

Solution x = [-0.19380928  0.36860987  0.63397973]
Residual = 2.7755575615628914e-17
In [52]:
# Question 1
A1 = np.array([[2, -1, 1],
               [1,  1, 1],
               [3, -1, -1]], float)
b1 = np.array([4, 3, 1], float)
x1 = mygauss(A1, b1)
print("\nQuestion 1: solution x =", x1)

Question 1: solution x = [1. 0. 2.]
In [53]:
# Question 2
A2 = np.array([[1.0, 5.0],
               [1.5, 7.501]], float)
b2 = np.array([17.0, 25.503], float)
x2 = mygauss(A2, b2)
print("\nQuestion 2: solution x =", x2)

Question 2: solution x = [2. 3.]

Question 4¶

Write your own function to compute the exponential function directly from the Taylor series:

$$ e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dots $$

The series should end when the last term is less than $10^{-6}$. Test your function against the built-in function exp, but be careful not to make $x$ too large — this could cause rounding error.

In [2]:
def calculate_exp(x):
    # 1e-6 is used as the stopping condition

    term = 1.0
    sum = term
    n = 1
    while abs(term) > 1e-6:
        term *= x/n
        sum += term
        n += 1
    return sum
In [4]:
test_values = [-1, 0, 1, 2, 5, 10, 20]
print("x\tdifference\tcalculate_exp(x)\t math.exp(x)")
for x in test_values:
    my_val = calculate_exp(x)
    true_val = math.exp(x)
    print(f"{x}\t{abs(my_val - true_val):.6e}\t{my_val:.6f}\t{true_val:.6f}")

x	difference	calculate_exp(x)	 math.exp(x)
-1	2.311427e-08	0.367879	0.367879
0	0.000000e+00	1.000000	1.000000
1	2.731266e-08	2.718282	2.718282
2	2.860474e-08	7.389056	7.389056
5	1.198070e-07	148.413159	148.413159
10	2.563938e-07	22026.465795	22026.465795
20	0.000000e+00	485165195.409790	485165195.409790

Question 5¶

If a random variable $X$ is distributed normally with zero mean and unit standard deviation, the probability that $0 \leq X \leq x$ is given by the standard normal function $\Phi(x)$. This is usually looked up in tables, but it may be approximated as follows:

$$ \Phi(x) = 0.5 - r\bigl(a t + b t^2 + c t^3\bigr), $$

where $$ a = 0.4361836,\quad b = -0.1201676,\quad c = 0.937298, $$ $$ r = \frac{\exp(-0.5 x^2)}{\sqrt{2\pi}},\quad t = \frac{1}{1 + 0.3326\,x}. $$

Write a function to compute $\Phi(x)$, and use it in a program to write out its values for $0 \leq x \leq 4$ in steps of 0.1. Check: $\Phi(1) = 0.3413.$

In [5]:
def calculate_phi_function(x):
    a = 0.4361836
    b = -0.1201676
    c = 0.937298
    r = np.exp(-0.5 * x**2) / np.sqrt(2 * np.pi)
    t = 1 / (1 + 0.3326 * x)
    return 0.5 - r * (a*t + b*t**2 + c*t**3)
In [6]:
x_values = np.arange(0, 4.1, 0.1)

print("x\tcalculate_phi_function(x)")
for x in x_values:
    print(f"{x:.1f}\t{calculate_phi_function(x):.4f}")

x	calculate_phi_function(x)
0.0	0.0000
0.1	0.0398
0.2	0.0793
0.3	0.1179
0.4	0.1554
0.5	0.1914
0.6	0.2257
0.7	0.2580
0.8	0.2881
0.9	0.3159
1.0	0.3413
1.1	0.3643
1.2	0.3849
1.3	0.4032
1.4	0.4192
1.5	0.4332
1.6	0.4452
1.7	0.4554
1.8	0.4641
1.9	0.4713
2.0	0.4772
2.1	0.4821
2.2	0.4861
2.3	0.4893
2.4	0.4918
2.5	0.4938
2.6	0.4953
2.7	0.4965
2.8	0.4974
2.9	0.4981
3.0	0.4986
3.1	0.4990
3.2	0.4993
3.3	0.4995
3.4	0.4997
3.5	0.4998
3.6	0.4998
3.7	0.4999
3.8	0.4999
3.9	0.5000
4.0	0.5000
In [8]:
print("\ncheck for x = 1: Φ(1) =", round(calculate_phi_function(1),4))

check for x = 1: Φ(1) = 0.3413

Question 6.a¶

There are many formulae for computing π (the ratio of a circle’s circumference to its diameter).
The simplest is:

$$ \frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \dots \tag{1} $$

which comes from putting $x = 1$ in the series:

$$ \arctan x = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \frac{x^9}{9} - \dots \tag{2} $$

a. Write a program to compute π using Equation (1). Use as many terms in the series as your computer will reasonably allow (start modestly, with 100 terms, say, and re-run your program with more and more each time). You should find that the series converges very slowly, i.e., it takes a lot of terms to get fairly close to π.

In [12]:
def calculate_pi_by_given_series(n_terms):
    terms = np.arange(n_terms)
    series = (-1)**terms / (2*terms + 1)
    return 4 * np.sum(series)
In [16]:
n_terms = 1000
pi_a = calculate_pi_by_given_series(n_terms)
print("6.a) calculated pi using given series: {:.10f}".format(pi_a))

6.a) calculated pi using given series: 3.1405926538

Question 6.b¶

b. Rearranging the series speeds up the convergence:

$$ \frac{\pi}{8} = \frac{1}{1 \times 3} + \frac{1}{5 \times 7} + \frac{1}{9 \times 11} + \dots $$

Write a program to compute π using this series instead. You should find that you need fewer terms to reach the same level of accuracy that you got in (a).

In [17]:
def calculate_pi_by_fast_series(n_terms):
    k = np.arange(n_terms)
    series = 1 / ((4*k + 1) * (4*k + 3))
    return 8 * np.sum(series)
In [18]:
n_terms_b = 100
pi_b = calculate_pi_by_fast_series(n_terms_b)
print("6.b) calculated pi using given faster series: {:.10f}".format(pi_b))

6.b) calculated pi using given faster series: 3.1365926848

Question 6.c¶

c. One of the fastest series for π is:

$$ \frac{\pi}{4} = 6 \arctan \frac{1}{8} + 2 \arctan \frac{1}{57} + \arctan \frac{1}{239} $$

Use this formula to compute π. Don’t use the built-in function atan to compute the arctangents, since that would be cheating. Rather use Equation (2).

In [19]:
def calculate_pi_by_arctan_series(x, n_terms=100):
    k = np.arange(n_terms)
    series = (-1)**k * x**(2*k + 1) / (2*k + 1)
    return np.sum(series)
In [21]:
pi_c = 4 * (6*calculate_pi_by_arctan_series(1/8) + 2*calculate_pi_by_arctan_series(1/57) + calculate_pi_by_arctan_series(1/239))
print("6.c) calculated pi using arctan series: {:.10f}".format(pi_c))

6.c) calculated pi using arctan series: 3.1415926536

Question 6.d¶

d. Can you vectorize any of your solutions (if you haven’t already)?

In [23]:
print("All the calculations above are vectorized by using numpy arrays.")

All the calculations above are vectorized by using numpy arrays.

Question 7¶

The following method of computing π is due to Archimedes:

  1. Let (A = 1) and (N = 6)
  2. Repeat 10 times, say:
    • Replace (N) by (2N)
    • Replace (A) by $$ A \gets \sqrt{2 - \sqrt{4 - A^2}} $$
    • Let $$ L = \frac{N \cdot A}{2} $$
    • Let $$ U = \frac{L}{\sqrt{1 - A^2/2}} $$
    • Let $$ P = \frac{U + L}{2} \quad \text{(estimate of } \pi\text{)} $$
    • Let $$ E = \frac{U - L}{2} \quad \text{(estimate of error)} $$
    • Print (N), (P), (E)
  3. Stop.

Write a program to implement the algorithm.

In [27]:
A = 1.0
N = 6

print("Iter\tN\tP (pi estimate)\tE (error estimate)")
for i in range(10):
    N = 2 * N
    
    A = math.sqrt(2 - math.sqrt(4 - A**2))

    L = N * A / 2
    U = L / (math.sqrt(1 - A**2 / 2))
    
    P = (U + L) / 2
    E = (U - L) / 2

    print(f"{i+1}\t{N}\t{P:.8f}\t{E:.8f}")

Iter	N	P (pi estimate)	E (error estimate)
1	12	3.22162925	0.11580071
2	24	3.16001597	0.02738736
3	48	3.14610799	0.00675779
4	96	3.14271595	0.00168400
5	192	3.14187313	0.00042066
6	384	3.14166275	0.00010514
7	768	3.14161018	0.00002628
8	1536	3.14159703	0.00000657
9	3072	3.14159375	0.00000164
10	6144	3.14159293	0.00000041

Question 8¶

If an amount of money $A$ is invested for $k$ years at a nominal annual interest rate $r$ (expressed as a decimal fraction), the value $V$ of the investment after $k$ years is given by:

$$ V = A \left(1 + \frac{r}{n}\right)^{n k} $$

where $n$ is the number of compounding periods per year.

Write a program to compute $V$ as $n$ gets larger and larger, i.e., as the compounding periods become more and more frequent, like monthly, daily, hourly, etc. Take $A = 1000$, $r = 0.04$ (4 per cent) and $k = 10$ years. You should observe that your output gradually approaches a limit.

Hint: Use a for loop which doubles (n) each time, starting with $n = 1$.

Also compute the value of the formula $A e^{r k}$ for the same values of $A$, $r$, and $k$ (use the exp function), and compare this value with the values of $V$ computed above. What do you conclude?

In [28]:
A = 1000
r = 0.04  
k = 10

n = 1
print("n\tV (compounded)\tV (continuous)")
for i in range(10):
    V = A * (1 + r/n)**(n*k)
    V_cont = A * np.exp(r*k)
    print(f"{n}\t{V:.6f}\t{V_cont:.6f}")
    n *= 2

print("\nContinuous compounding value: {:.6f}".format(V_cont))

n	V (compounded)	V (continuous)
1	1480.244285	1491.824698
2	1485.947396	1491.824698
4	1488.863734	1491.824698
8	1490.338568	1491.824698
16	1491.080212	1491.824698
32	1491.452099	1491.824698
64	1491.638309	1491.824698
128	1491.731481	1491.824698
256	1491.778084	1491.824698
512	1491.801389	1491.824698

Continuous compounding value: 1491.824698

  • As the number of compounding periods n increases, the value V of the investment converges to the calculated continuous compounding value.

  • This shows that the continous compunding value is an upper limit for compounded value and this value approaches to the continous compounding value as the compounding frequence increases.