Numerical Analysis

Numerical Analysis

This module collects numerical methods I implemented in C++ and exposed through Python.

Source Code

The source code is available on GitHub.

Package Documentation

The available functions and examples are documented below.

Main Reference

The main reference is Numerical Analysis, 9th edition, by Richard L. Burden and J. Douglas Faires.

---

Why Use Numerical Analysis?

Numerical analysis uses approximation algorithms to solve mathematical problems when a symbolic solution is unavailable or impractical.

Example Applications

• Example 1

  Suppose we manage a city whose population at time $t$ is $N(t)$, with net birth rate $\lambda$ and immigration rate $v$. We can model the population with the differential equation

  $$\dfrac{dN(t)}{dt} = \lambda N(t) + v$$

  Its solution is

  $$N(t) = N(0)e^{\lambda t} + \dfrac{v}{\lambda} (e ^ {\lambda t} - 1)$$

  Suppose the city had a population of 10,000 last year and received 2,000 immigrants. Its current population is 15,500, giving us

  $$15,500 = 10,000e^{\lambda} + \dfrac{2,000}{\lambda} (e ^ {\lambda} - 1)$$

  We need to find the net birth rate $\lambda$, but we cannot isolate it algebraically. A numerical root-finding method gives us an approximation.

• Example 2

  Suppose a store's sales depend on its local advertising spend. We can write

  $$Sales = f(Advertising)$$

  If the store spends $50,000 on advertising, we may want to estimate how sales change as spending changes:

  $$\dfrac{dSales}{dAdvertising} = f'(Advertising); \enspace Advertising = 50,000$$

  Without the functional form of $f(Advertising)$, we must approximate the derivative from observed data.

---

Documentation

• Solutions of Single Variable Equations

  Approximate $x$ such that $f(x) = y, \enspace x, y \in \R$

  Implementation:

  The implementation uses the secant method, which generally converges faster than bisection when its assumptions hold.

  Implementation Example

  Solve for $x$ in

  1. $14 = x^2 + 2x - 1$
  2. $0 = x^2 + 1$

  Exact Solution:

  1. $x = 3$
  2. $x \notin \R$

  Approximation:

  >>> from pysiclib import numerical
  >>> def example_func(x):
  ... 	return x ** 2 + 2 * x - 1
  >>> target_y = 14
  >>> numerical.equation_solution(example_func, target_y)
  3.0
  >>> def no_solution_example(x):
  ... 	return x ** 2 + 1
  >>> target_y = 0
  >>> numerical.equation_solution(no_solution_example, target_y)
  None

• Differentiation

  Approximate $\dfrac{d f(x)}{d x}$ given an array $X = [(x_1, f(x_1)), ..., (x_n, f(x_n))]$.

  Implementation:

  The implementation uses the five-point midpoint method when $(x, f(x)) \in X$.

  Implementation Example

  Solve $f'(5), \enspace f(x) = x^2 + x$

  Exact Solution: $f'(5) = 11$

  Approximation:

  >>> from pysiclib import numerical, linalg
  >>> data = linalg.Tensor([x ** 2 + x for x in range(11)])
  >>> x_index = 5
  >>> numerical.derivative_at_index(data, x_index)
  11.0

• Integration

  Approximate $F(x),\enspace F(x) = \int_{a}^{b}f(x)dx$ given an array $X = [(x_1, f(x_1)),..., (x_n, f(x_n))]$.

  Implementation:

  The implementation uses composite Simpson's rule.

  Implementation Example

  Solve $I = \int_{0}^{5}(x^2 + x)dx$

  Exact Solution: $I = 325/6 = 51.1666...$

  Approximation:

  >>> from pysiclib import numerical, linalg
  >>> unit_steps = 100
  >>> data = linalg.Tensor(
  ...			[(x / unit_steps) ** 2 + (x / unit_steps) for x in
  ...			range(11 * unit_steps)])]
  ...		)
  >>> int_start, int_end = 0, 5 * unit_steps
  ... #note the interval [start, end] is integrated over
  >>> numerical.integral_index_interval(data, int_start, int_end)
  54.107366

• Differential Equations - Initial Value problems

  For given $t \in [a, b]$ approximate $y_1(t), ..., y_n(t)$ given an $n$-order system of initial value problems having the form

  $$\dfrac{dy_1}{dt} = f_1(t, y_1, ..., y_n)\\ \dfrac{dy_2}{dt} = f_2(t, y_1, ..., y_n)\\ .\\ .\\ .\\ \dfrac{dy_n}{dt} = f_n(t, y_1, ..., y_n)\\$$

  With an initial condition of $\lambda \in [a, b]$ such that

  $$y_1(\lambda) = \alpha_1\\ y_2(\lambda) = \alpha_2\\ .\\ .\\ .\\ y_n(\lambda) = \alpha_n$$

  Where there may exist an arbitrary amount of $y_i(t) = y'_j(t),\enspace i \neq j$

  Implementation

  The implementation uses a Runge-Kutta method.

  Implementation Example

  $$\begin{aligned} y_1' &= -4y_1 + 3y_2 + 6, \enspace &y_1(0) = 0\\ y_2' &= 0.6y'_1 - 0.2y_2, \enspace &y_2(0) = 0 \end{aligned}$$

  Let $Y(t) = (y_1(t), y_2(t))$, find $Y(0.5)$

  Exact Solution: $Y(0.5) = (1.793527048, 1.014415452)$

  Approximation:

  >>> from pysiclib import numerical, linalg
  >>> def system_of_eqs(t, var_arr):
  ...   var_arr = var_arr.get_buffer()[:] #temporary for now
  ... 	dvar_arr = [0.0, 0.0]
  ... 	dvar_arr[0] = -4 * var_arr[0] + 3 * var_arr[1] + 6
  ... 	dvar_arr[1] = 0.6 * dvar_arr[0] - 0.2 * var_arr[1]
  ... 	return linalg.Tensor(dvar_arr)
  >>> init_cond = linalg.Tensor([0.0, 0.0])
  >>> init_val = 0.0
  >>> target_val = 0.5
  >>> numerical.initial_value_problem(
  ... 	system_of_eqs, init_cond, target_val, init_val)
  [1.79352705 1.01441545]