Numerical Analysis

01

Bisection Method

f(x) = x³ - x - 2 = 0, interval [1.0, 2.0]

expand_more

code 01_bisection.py

def bisection(f, a, b, tol=1e-5, max_iter=50):
    if f(a) * f(b) >= 0:
        return print("Error: f(a) and f(b) must have opposite signs.")

    print(f"{'Iter':<6}{'a':<12}{'b':<12}{'c (Mid)':<14}{'f(c)':<12}")
    print("-" * 56)

    for i in range(1, max_iter + 1):
        c = (a + b) / 2
        fc = f(c)
        print(f"{i:<6}{a:<12.5f}{b:<12.5f}{c:<14.5f}{fc:<12.5f}")

        if abs(fc) < tol or (b - a) / 2 < tol:
            return c
        a, b = (a, c) if f(a) * fc < 0 else (c, b)

    return (a + b) / 2

def f(x):
    return x**3 - x - 2

root = bisection(f, 1.0, 2.0)
print(f"\nRoot: {root:.5f}")

terminal Terminal Output Trace

python3 01_bisection.py

Iter a b c (Mid) f(c) -------------------------------------------------------- 1 1.00000 2.00000 1.50000 -0.12500 2 1.50000 2.00000 1.75000 1.60938 3 1.50000 1.75000 1.62500 0.66602 4 1.50000 1.62500 1.56250 0.25220 5 1.50000 1.56250 1.53125 0.05911 6 1.50000 1.53125 1.51562 -0.03405 7 1.51562 1.53125 1.52344 0.01225 8 1.51562 1.52344 1.51953 -0.01097 9 1.51953 1.52344 1.52148 0.00062 10 1.51953 1.52148 1.52051 -0.00518 11 1.52051 1.52148 1.52100 -0.00228 12 1.52100 1.52148 1.52124 -0.00083 13 1.52124 1.52148 1.52136 -0.00010 14 1.52136 1.52148 1.52142 0.00026 15 1.52136 1.52142 1.52139 0.00008 16 1.52136 1.52139 1.52138 -0.00001 17 1.52138 1.52139 1.52139 0.00003 Root: 1.52139

02

Regula Falsi (False Position)

f(x) = x³ - 2x - 5 = 0, interval [2.0, 3.0]

expand_more

code 02_regula_falsi.py

def regula_falsi(f, a, b, tol=1e-5, max_iter=50):
    if f(a) * f(b) >= 0:
        return print("Error: f(a) and f(b) must have opposite signs.")

    print(f"{'Iter':<6}{'a':<12}{'b':<12}{'c':<14}{'f(c)':<12}")
    print("-" * 56)

    for i in range(1, max_iter + 1):
        c = (a * f(b) - b * f(a)) / (f(b) - f(a))
        fc = f(c)
        print(f"{i:<6}{a:<12.5f}{b:<12.5f}{c:<14.5f}{fc:<12.5f}")

        if abs(fc) < tol:
            return c
        a, b = (a, c) if f(a) * fc < 0 else (c, b)

    return c

def f(x):
    return x**3 - 2 * x - 5

root = regula_falsi(f, 2.0, 3.0)
print(f"\nRoot: {root:.5f}")

terminal Terminal Output Trace

python3 02_regula_falsi.py

Iter a b c f(c) -------------------------------------------------------- 1 2.00000 3.00000 2.05882 -0.39080 2 2.05882 3.00000 2.08126 -0.14720 3 2.08126 3.00000 2.08964 -0.05468 4 2.08964 3.00000 2.09274 -0.02020 5 2.09274 3.00000 2.09388 -0.00745 6 2.09388 3.00000 2.09431 -0.00275 7 2.09431 3.00000 2.09446 -0.00101 8 2.09446 3.00000 2.09452 -0.00037 9 2.09452 3.00000 2.09454 -0.00014 10 2.09454 3.00000 2.09455 -0.00005 11 2.09455 3.00000 2.09455 -0.00002 12 2.09455 3.00000 2.09455 -0.00001 Root: 2.09455

03

Newton Raphson Method

f(x) = x³ - 3x + 1 = 0, initial guess x₀ = 1.5

expand_more

code 03_newton_raphson.py

def newton_raphson(f, df, x0, tol=1e-5, max_iter=50):
    print(f"{'Iter':<6}{'x0':<12}{'f(x0)':<14}{'df(x0)':<14}{'x1':<12}")
    print("-" * 58)

    for i in range(1, max_iter + 1):
        fx, dfx = f(x0), df(x0)
        if dfx == 0:
            return print("Derivative is zero.")
        x1 = x0 - fx / dfx
        print(f"{i:<6}{x0:<12.5f}{fx:<14.5f}{dfx:<14.5f}{x1:<12.5f}")

        if abs(x1 - x0) < tol or abs(f(x1)) < tol:
            return x1
        x0 = x1

    return x0

def f(x):
    return x**3 - 3 * x + 1

def df(x):
    return 3 * x**2 - 3

root = newton_raphson(f, df, 1.5)
print(f"\nRoot: {root:.5f}")

terminal Terminal Output Trace

python3 03_newton_raphson.py

Iter x0 f(x0) df(x0) x1 ---------------------------------------------------------- 1 1.50000 -0.12500 3.75000 1.53333 2 1.53333 0.00504 4.05333 1.53209 Root: 1.53209

04

Secant Method

f(x) = x² - 4 = 0, x₀ = 1.0, x₁ = 3.0

expand_more

code 04_secant.py

def secant(f, x0, x1, tol=1e-5, max_iter=50):
    print(f"{'Iter':<6}{'x0':<12}{'x1':<12}{'f(x1)':<14}{'x2':<12}")
    print("-" * 56)

    for i in range(1, max_iter + 1):
        f0, f1 = f(x0), f(x1)
        if f1 - f0 == 0:
            break
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        print(f"{i:<6}{x0:<12.5f}{x1:<12.5f}{f1:<14.5f}{x2:<12.5f}")

        if abs(x2 - x1) < tol or abs(f(x2)) < tol:
            return x2
        x0, x1 = x1, x2

    return x1

def f(x):
    return x**2 - 4

root = secant(f, 1.0, 3.0)
print(f"\nRoot: {root:.5f}")

terminal Terminal Output Trace

python3 04_secant.py

Iter x0 x1 f(x1) x2 -------------------------------------------------------- 1 1.00000 3.00000 5.00000 1.75000 2 3.00000 1.75000 -0.93750 1.94737 3 1.75000 1.94737 -0.20776 2.00356 4 1.94737 2.00356 0.01425 1.99995 5 2.00356 1.99995 -0.00019 2.00000 Root: 2.00000

05

Fixed Point Iteration

g(x) = cos(x), x₀ = 0.5

expand_more

code 05_fixed_point.py

import math

def fixed_point(g, x0, tol=1e-5, max_iter=50):
    print(f"{'Iter':<6}{'x0':<14}{'x1 = g(x0)':<14}")
    print("-" * 34)

    for i in range(1, max_iter + 1):
        x1 = g(x0)
        print(f"{i:<6}{x0:<14.5f}{x1:<14.5f}")
        if abs(x1 - x0) < tol:
            return x1
        x0 = x1

    return x0

def g(x):
    return math.cos(x)

root = fixed_point(g, 0.5)
print(f"\nFixed Point: {root:.5f}")

terminal Terminal Output Trace

python3 05_fixed_point.py

Iter x0 x1 = g(x0) ---------------------------------- 1 0.50000 0.87758 2 0.87758 0.63901 3 0.63901 0.80269 4 0.80269 0.69478 5 0.69478 0.76820 6 0.76820 0.71917 7 0.71917 0.75236 8 0.75236 0.73008 9 0.73008 0.74512 10 0.74512 0.73501 11 0.73501 0.74183 12 0.74183 0.73724 13 0.73724 0.74033 14 0.74033 0.73825 15 0.73825 0.73965 16 0.73965 0.73870 17 0.73870 0.73934 18 0.73934 0.73891 19 0.73891 0.73920 20 0.73920 0.73901 21 0.73901 0.73914 22 0.73914 0.73905 23 0.73905 0.73911 24 0.73911 0.73907 25 0.73907 0.73910 26 0.73910 0.73908 27 0.73908 0.73909 28 0.73909 0.73908 Fixed Point: 0.73908

06

Gauss Elimination Method

Forward elimination and back substitution

expand_more

code 06_gauss_elimination.py

def gauss_elimination(A, b):
    n = len(b)
    # Forward Elimination
    for i in range(n):
        for j in range(i + 1, n):
            ratio = A[j][i] / A[i][i]
            for k in range(i, n):
                A[j][k] -= ratio * A[i][k]
            b[j] -= ratio * b[i]

    # Back Substitution
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - sum(A[i][j] * x[j] for j in range(i + 1, n))) / A[i][i]
    return x

A = [[2.0, 1.0, -1.0], [-3.0, -1.0, 2.0], [-2.0, 1.0, 2.0]]
b = [8.0, -11.0, -3.0]
ans = gauss_elimination(A, b)
for i, val in enumerate(ans):
    print(f"x{i+1} = {val:.4f}")

terminal Terminal Output Trace

python3 06_gauss_elimination.py

x1 = 2.0000 x2 = 3.0000 x3 = -1.0000

07

Gauss-Jordan Method

Diagonalization to reduced row echelon form

expand_more

code 07_gauss_jordan.py

def gauss_jordan(A, b):
    n = len(b)
    for i in range(n):
        p = A[i][i]
        A[i] = [v / p for v in A[i]]
        b[i] /= p
        for k in range(n):
            if k != i:
                factor = A[k][i]
                A[k] = [A[k][j] - factor * A[i][j] for j in range(n)]
                b[k] -= factor * b[i]
    return b

A = [[2.0, 3.0, 1.0], [1.0, 2.0, 3.0], [3.0, 1.0, 2.0]]
b = [9.0, 6.0, 8.0]
ans = gauss_jordan(A, b)
for i, val in enumerate(ans):
    print(f"x{i+1} = {val:.4f}")

terminal Terminal Output Trace

python3 07_gauss_jordan.py

x1 = 1.9444 x2 = 1.6111 x3 = 0.2778

08

Gauss-Jacobi Iteration

Simultaneous displacement iterative solver

expand_more

code 08_gauss_jacobi.py

def gauss_jacobi(A, b, tol=1e-5, max_iter=50):
    n = len(b)
    x = [0.0] * n
    print(f"{'Iter':<6}" + "".join(f"{'x' + str(i+1):<12}" for i in range(n)))
    print("-" * (6 + 12 * n))

    for step in range(1, max_iter + 1):
        x_new = [(b[i] - sum(A[i][j] * x[j] for j in range(n) if j != i)) / A[i][i] for i in range(n)]
        print(f"{step:<6}" + "".join(f"{v:<12.5f}" for v in x_new))

        if max(abs(x_new[i] - x[i]) for i in range(n)) < tol:
            return x_new
        x = x_new
    return x

A = [[10.0, -1.0, 2.0], [-1.0, 11.0, -1.0], [2.0, -1.0, 10.0]]
b = [6.0, 25.0, -11.0]
sol = gauss_jacobi(A, b)
print("\nSolution:", [round(v, 4) for v in sol])

terminal Terminal Output Trace

python3 08_gauss_jacobi.py

Iter x1 x2 x3 ------------------------------------------ 1 0.60000 2.27273 -1.10000 2 1.04727 2.22727 -0.99273 3 1.02127 2.27769 -1.08673 4 1.04511 2.26678 -1.07649 5 1.04197 2.26988 -1.08235 6 1.04346 2.26906 -1.08141 7 1.04319 2.26928 -1.08179 8 1.04328 2.26922 -1.08171 9 1.04326 2.26923 -1.08174 10 1.04327 2.26923 -1.08173 Solution: [1.0433, 2.2692, -1.0817]

09

Gauss-Seidel Iteration

Successive displacement in-place solver

expand_more

code 09_gauss_seidel.py

def gauss_seidel(A, b, tol=1e-5, max_iter=50):
    n = len(b)
    x = [0.0] * n
    print(f"{'Iter':<6}" + "".join(f"{'x' + str(i+1):<12}" for i in range(n)))
    print("-" * (6 + 12 * n))

    for step in range(1, max_iter + 1):
        x_old = list(x)
        for i in range(n):
            s = sum(A[i][j] * x[j] for j in range(n) if j != i)
            x[i] = (b[i] - s) / A[i][i]

        print(f"{step:<6}" + "".join(f"{v:<12.5f}" for v in x))
        if max(abs(x[i] - x_old[i]) for i in range(n)) < tol:
            return x
    return x

A = [[10.0, -1.0, 2.0], [-1.0, 11.0, -1.0], [2.0, -1.0, 10.0]]
b = [6.0, 25.0, -11.0]
sol = gauss_seidel(A, b)
print("\nSolution:", [round(v, 4) for v in sol])

terminal Terminal Output Trace

python3 09_gauss_seidel.py

Iter x1 x2 x3 ------------------------------------------ 1 0.60000 2.32727 -0.98727 2 1.03018 2.27663 -1.07837 3 1.04334 2.26954 -1.08171 4 1.04330 2.26923 -1.08174 5 1.04327 2.26923 -1.08173 6 1.04327 2.26923 -1.08173 Solution: [1.0433, 2.2692, -1.0817]

10

Lagrange Interpolation

Arbitrary/unequal nodal spacing formula

expand_more

code 10_lagrange.py

import math

def lagrange(x_vals, y_vals, target_x):
    n = len(x_vals)
    print(f"{'i':<6}{'Basis L_i(x)':<16}{'y_i':<10}{'Term':<12}")
    print("-" * 44)

    total = 0.0
    for i in range(n):
        L_i = math.prod((target_x - x_vals[j]) / (x_vals[i] - x_vals[j]) for j in range(n) if i != j)
        term = y_vals[i] * L_i
        total += term
        print(f"{i:<6}{L_i:<16.5f}{y_vals[i]:<10.2f}{term:<12.5f}")

    return total

ans = lagrange([5.0, 6.0, 9.0, 11.0], [12.0, 13.0, 14.0, 16.0], 10.0)
print(f"\nInterpolated value: {ans:.5f}")

terminal Terminal Output Trace

python3 10_lagrange.py

i Basis L_i(x) y_i Term -------------------------------------------- 0 0.16667 12.00 2.00000 1 -0.33333 13.00 -4.33333 2 0.83333 14.00 11.66667 3 0.33333 16.00 5.33333 Interpolated value: 14.66667

11

Newton's Divided Difference

Divided difference table interpolation

expand_more

code 11_newton_divided_diff.py

def newton_divided_diff(x, y, target_x):
    n = len(x)
    table = [[0.0] * n for _ in range(n)]
    for i in range(n):
        table[i][0] = y[i]

    for j in range(1, n):
        for i in range(n - j):
            table[i][j] = (table[i + 1][j - 1] - table[i][j - 1]) / (x[i + j] - x[i])

    print("Divided Difference Table:")
    for i in range(n):
        print(f"{x[i]:<6.1f}" + "".join(f"{table[i][j]:<12.4f}" for j in range(n - i)))

    res, prod_term = table[0][0], 1.0
    for j in range(1, n):
        prod_term *= (target_x - x[j - 1])
        res += table[0][j] * prod_term
    return res

ans = newton_divided_diff([4.0, 5.0, 7.0, 10.0, 11.0], [48.0, 100.0, 294.0, 900.0, 1210.0], 6.0)
print(f"\nInterpolated Value: {ans:.5f}")

terminal Terminal Output Trace

python3 11_newton_divided_diff.py

Divided Difference Table: 4.0 48.0000 52.0000 15.0000 1.0000 0.0000 5.0 100.0000 97.0000 21.0000 1.0000 7.0 294.0000 202.0000 27.0000 10.0 900.0000 310.0000 11.0 1210.0000 Interpolated Value: 180.00000

12

Newton's Forward Interpolation

Forward finite difference formula

expand_more

code 12_newton_forward.py

import math

def newton_forward(x, y, target_x):
    n = len(x)
    h, u = x[1] - x[0], (target_x - x[0]) / (x[1] - x[0])
    diff = [[0.0] * n for _ in range(n)]
    for i in range(n):
        diff[0][i] = y[i]

    for order in range(1, n):
        for i in range(n - order):
            diff[order][i] = diff[order - 1][i + 1] - diff[order - 1][i]

    print("Forward Difference Table:")
    for i in range(n):
        print(f"{x[i]:<6.1f}" + "".join(f"{diff[j][i]:<12.4f}" for j in range(n - i)))

    res, u_term = diff[0][0], 1.0
    for order in range(1, n):
        u_term *= (u - (order - 1))
        res += (u_term * diff[order][0]) / math.factorial(order)
    return res

ans = newton_forward([10.0, 20.0, 30.0, 40.0, 50.0], [0.1736, 0.3420, 0.5000, 0.6428, 0.7660], 15.0)
print(f"\nInterpolated Value: {ans:.5f}")

terminal Terminal Output Trace

python3 12_newton_forward.py

Forward Difference Table: 10.0 0.1736 0.1684 -0.0104 -0.0048 0.0004 20.0 0.3420 0.1580 -0.0152 -0.0044 30.0 0.5000 0.1428 -0.0196 40.0 0.6428 0.1232 50.0 0.7660 Interpolated Value: 0.25878

13

Trapezoidal Rule

∫₀¹ 1/(1 + x²) dx with n = 6 intervals

expand_more

code 13_trapezoidal.py

def trapezoidal(f, a, b, n=6):
    h = (b - a) / n
    x_pts = [a + i * h for i in range(n + 1)]
    y_pts = [f(x) for x in x_pts]

    print(f"{'i':<6}{'x_i':<12}{'f(x_i)':<14}")
    print("-" * 32)
    for i, (xi, yi) in enumerate(zip(x_pts, y_pts)):
        print(f"{i:<6}{xi:<12.4f}{yi:<14.5f}")

    integral = (h / 2.0) * (y_pts[0] + 2 * sum(y_pts[1:n]) + y_pts[n])
    return integral

def f(x):
    return 1 / (1 + x**2)

ans = trapezoidal(f, 0.0, 1.0, n=6)
print(f"\nIntegral: {ans:.5f}")

terminal Terminal Output Trace

python3 13_trapezoidal.py

i x_i f(x_i) -------------------------------- 0 0.0000 1.00000 1 0.1667 0.97297 2 0.3333 0.90000 3 0.5000 0.80000 4 0.6667 0.69231 5 0.8333 0.59016 6 1.0000 0.50000 Integral: 0.78424

14

Simpson's 1/3 Rule

∫₀¹ 1/(1 + x²) dx with n = 6 intervals

expand_more

code 14_simpson_1_3.py

def simpson_1_3(f, a, b, n=6):
    if n % 2 != 0:
        n += 1
    h = (b - a) / n
    x_pts = [a + i * h for i in range(n + 1)]
    y_pts = [f(x) for x in x_pts]

    print(f"{'i':<6}{'x_i':<12}{'f(x_i)':<14}{'Mult':<6}")
    print("-" * 38)
    for i, (xi, yi) in enumerate(zip(x_pts, y_pts)):
        mult = 1 if i in (0, n) else (4 if i % 2 != 0 else 2)
        print(f"{i:<6}{xi:<12.4f}{yi:<14.5f}{mult:<6}")

    integral = (h / 3.0) * (y_pts[0] + y_pts[n] + 4 * sum(y_pts[1:n:2]) + 2 * sum(y_pts[2:n:2]))
    return integral

def f(x):
    return 1 / (1 + x**2)

ans = simpson_1_3(f, 0.0, 1.0, n=6)
print(f"\nIntegral: {ans:.5f}")

terminal Terminal Output Trace

python3 14_simpson_1_3.py

i x_i f(x_i) Mult -------------------------------------- 0 0.0000 1.00000 1 1 0.1667 0.97297 4 2 0.3333 0.90000 2 3 0.5000 0.80000 4 4 0.6667 0.69231 2 5 0.8333 0.59016 4 6 1.0000 0.50000 1 Integral: 0.78540

15

Simpson's 3/8 Rule

∫₀¹ 1/(1 + x²) dx with n = 6 intervals

expand_more

code 15_simpson_3_8.py

def simpson_3_8(f, a, b, n=6):
    while n % 3 != 0:
        n += 1
    h = (b - a) / n
    x_pts = [a + i * h for i in range(n + 1)]
    y_pts = [f(x) for x in x_pts]

    print(f"{'i':<6}{'x_i':<12}{'f(x_i)':<14}{'Mult':<6}")
    print("-" * 38)
    for i, (xi, yi) in enumerate(zip(x_pts, y_pts)):
        mult = 1 if i in (0, n) else (2 if i % 3 == 0 else 3)
        print(f"{i:<6}{xi:<12.4f}{yi:<14.5f}{mult:<6}")

    middle = sum((2 if i % 3 == 0 else 3) * y_pts[i] for i in range(1, n))
    return (3 * h / 8.0) * (y_pts[0] + y_pts[n] + middle)

def f(x):
    return 1 / (1 + x**2)

ans = simpson_3_8(f, 0.0, 1.0, n=6)
print(f"\nIntegral: {ans:.5f}")

terminal Terminal Output Trace

python3 15_simpson_3_8.py

i x_i f(x_i) Mult -------------------------------------- 0 0.0000 1.00000 1 1 0.1667 0.97297 3 2 0.3333 0.90000 3 3 0.5000 0.80000 2 4 0.6667 0.69231 3 5 0.8333 0.59016 3 6 1.0000 0.50000 1 Integral: 0.78540

16

Euler's Method

dy/dx = x + y, y(0) = 1.0 (h = 0.05)

expand_more

code 16_euler.py

def euler(f, x0, y0, x_end, h=0.05):
    x, y, step = x0, y0, 0
    print(f"{'Step':<6}{'x':<10}{'y':<12}{'Slope':<12}")
    print("-" * 40)

    while x < x_end - 1e-9:
        slope = f(x, y)
        print(f"{step:<6}{x:<10.3f}{y:<12.5f}{slope:<12.5f}")
        y += h * slope
        x += h
        step += 1

    print(f"{step:<6}{x:<10.3f}{y:<12.5f}")
    return y

def f(x, y):
    return x + y

ans = euler(f, 0.0, 1.0, 0.2, h=0.05)
print(f"\ny(0.2): {ans:.5f}")

terminal Terminal Output Trace

python3 16_euler.py

Step x y Slope ---------------------------------------- 0 0.000 1.00000 1.00000 1 0.050 1.05000 1.10000 2 0.100 1.10500 1.20500 3 0.150 1.16525 1.31525 4 0.200 1.23101 y(0.2): 1.23101

17

Modified Euler's Method (Heun's)

Predictor-corrector scheme (h = 0.05)

expand_more

code 17_modified_euler.py

def modified_euler(f, x0, y0, x_end, h=0.05):
    x, y, step = x0, y0, 0
    print(f"{'Step':<6}{'x':<10}{'y':<12}{'y_predict':<14}")
    print("-" * 42)

    while x < x_end - 1e-9:
        k1 = f(x, y)
        y_predict = y + h * k1
        print(f"{step:<6}{x:<10.3f}{y:<12.5f}{y_predict:<14.5f}")
        y += (h / 2.0) * (k1 + f(x + h, y_predict))
        x += h
        step += 1

    print(f"{step:<6}{x:<10.3f}{y:<12.5f}")
    return y

def f(x, y):
    return x + y

ans = modified_euler(f, 0.0, 1.0, 0.2, h=0.05)
print(f"\ny(0.2): {ans:.5f}")

terminal Terminal Output Trace

python3 17_modified_euler.py

Step x y y_predict ------------------------------------------ 0 0.000 1.00000 1.05000 1 0.050 1.05250 1.10763 2 0.100 1.11025 1.17077 3 0.150 1.17353 1.23971 4 0.200 1.24261 y(0.2): 1.24261

18

Runge-Kutta 4th Order (RK4)

4th-order weighted slope predictor (h = 0.05)

expand_more

code 18_rk4.py

def rk4(f, x0, y0, x_end, h=0.05):
    x, y, step = x0, y0, 0
    print(f"{'Step':<6}{'x':<8}{'y':<10}{'k1':<9}{'k2':<9}{'k3':<9}{'k4':<9}")
    print("-" * 60)

    while x < x_end - 1e-9:
        k1 = h * f(x, y)
        k2 = h * f(x + h / 2, y + k1 / 2)
        k3 = h * f(x + h / 2, y + k2 / 2)
        k4 = h * f(x + h, y + k3)
        print(f"{step:<6}{x:<8.2f}{y:<10.4f}{k1:<9.4f}{k2:<9.4f}{k3:<9.4f}{k4:<9.4f}")
        y += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0
        x += h
        step += 1

    print(f"{step:<6}{x:<8.2f}{y:<10.4f}")
    return y

def f(x, y):
    return x + y

ans = rk4(f, 0.0, 1.0, 0.2, h=0.05)
print(f"\ny(0.2): {ans:.5f}")

terminal Terminal Output Trace

python3 18_rk4.py

Step x y k1 k2 k3 k4 ------------------------------------------------------------ 0 0.00 1.0000 0.0500 0.0525 0.0526 0.0551 1 0.05 1.0525 0.0551 0.0578 0.0578 0.0605 2 0.10 1.1103 0.0605 0.0633 0.0633 0.0662 3 0.15 1.1737 0.0662 0.0691 0.0692 0.0721 4 0.20 1.2428 y(0.2): 1.24281

19

Power Method

Dominant eigenvalue and eigenvector of A

expand_more

code 19_power_method.py

def power_method(A, tol=1e-5, max_iter=50):
    n, x, lambda_old = len(A), [1.0] * len(A), 0.0
    print(f"{'Iter':<6}{'Eigenvalue':<14}{'Eigenvector':<20}")
    print("-" * 40)

    for step in range(1, max_iter + 1):
        y = [sum(A[i][j] * x[j] for j in range(n)) for i in range(n)]
        lambda_new = max(y, key=abs)
        x = [v / lambda_new for v in y]
        vec_str = "[" + ", ".join(f"{v:.3f}" for v in x) + "]"
        print(f"{step:<6}{lambda_new:<14.4f}{vec_str:<20}")

        if abs(lambda_new - lambda_old) < tol:
            return lambda_new, x
        lambda_old = lambda_new

    return lambda_new, x

val, vec = power_method([[4.0, 1.0], [2.0, 3.0]])
print(f"\nEigenvalue: {val:.4f}, Eigenvector: {[round(v, 4) for v in vec]}")

terminal Terminal Output Trace

python3 19_power_method.py

Iter Eigenvalue Eigenvector ---------------------------------------- 1 5.0000 [1.000, 1.000] 2 5.0000 [1.000, 1.000] Eigenvalue: 5.0000, Eigenvector: [1.0, 1.0]

20

Linear Curve Fitting (Least Squares)

Best fit line y = mx + c from coordinate pairs

expand_more

code 20_linear_fit.py

def linear_fit(x, y):
    n = len(x)
    sx, sy = sum(x), sum(y)
    sxy = sum(a * b for a, b in zip(x, y))
    sx2 = sum(a**2 for a in x)

    print(f"n = {n}, sum_x = {sx:.2f}, sum_y = {sy:.2f}, sum_xy = {sxy:.2f}, sum_x2 = {sx2:.2f}")
    m = (n * sxy - sx * sy) / (n * sx2 - sx**2)
    c = (sy - m * sx) / n
    return m, c

m, c = linear_fit([1.0, 2.0, 3.0, 4.0, 5.0], [2.2, 2.8, 3.6, 4.5, 5.1])
print(f"Fitted Line: y = {m:.4f}x + {c:.4f}")

terminal Terminal Output Trace

python3 20_linear_fit.py

n = 5, sum_x = 15.00, sum_y = 18.20, sum_xy = 62.10, sum_x2 = 55.00 Fitted Line: y = 0.7500x + 1.3900