ECE 204: Numerical Methods¶

Linear regression¶

Given a regression $y=mx+b$ and a data set $(x_{i..n}, y_{i..n})$ , the residual is the difference between the actual and regressed data:

$E_i=y_i-b-mx_i$

Method of least squares¶

This method minimises the sum of the square of residuals.

$\boxed{S_r=\sum^n_{i=1}E_i^2}$

$m$ and $b$ can be found by taking the partial derivative and solving for them:

$\frac{\partial S_r}{\partial m}=0, \frac{\partial S_r}{\partial b}=0$

This returns, where $\overline y$ is the mean of the actual $y$ -values:

$\boxed{m=\frac{n\sum^n_{i=1}x_iy_i-\sum^n_{i=1}x_i\sum^n_{i=1}y_i}{n\sum^n_{i=1}x_i^2-\left(\sum^n_{i=1}x_i\right)^2}} \\ b=\overline y-m\overline x$

The total sum of square around the mean is based off of the actual data:

$\boxed{S_t=\sum(y_i-\overline y)^2}$

Error is measured with the coefficient of determination $r^2$ — the closer the value is to 1, the lower the error.

$r^2=\frac{S_t-S_r}{S_t}$

If the intercept is the origin, $m$ reduces down to a simpler form:

$m=\frac{\sum^n_{i=1}x_iy_i}{\sum^n_{i=1}x_i^2}$

Non-linear regression¶

Exponential regression¶

Solving for the same partial derivatives returns the same values, although bisection may be required for the exponent coefficient ( $e^{bx}$ ) Instead, linearising may make things easier (by taking the natural logarithm of both sides. Afterward, solving as if it were in the form $y=mx+b$ returns correct

Example

$y=ax^b\implies\ln y = \ln a + b\ln x$

Polynomial regression¶

The residiual is the offset at the end of a polynomial.

$y=a+bx+cx^2+E$

Taking the relevant partial derivatives returns a system of equations which can be solved in a matrix.

Interpolation¶

Interpolation ensures that every point is crossed.

Direct method¶

To interpolate $n+1$ data points, you need a polynomial of a degree up to $n$ , and points that enclose the desired value. Substituting the $x$ and $y$ values forms a system of equations for a polynomial of a degree equal to the number of points chosen - 1.

Newton's divided difference method¶

This method guesses the slope to interpolate. Where $x_0$ is an existing point:

$\boxed{f(x)=b_0+b_1(x-x_0)}$

The constant is an existing y-value and the slope is an average.

$\begin{align*} b_0&=f(x_0) \\ b_1&=\frac{f(x_1)-f(x_0)}{x_1-x_0} \end{align*}$

This extends to a quadratic, where the second slope is the average of the first two slopes:

$\boxed{f(x)=b_0+b_1(x-x_0)+b_2(x-x_0)(x-x_1)}$

$b_2=\frac{\frac{f(x_2)-f(x_1)}{x_2-x_1}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_0}$

Derivatives¶

Derivatives are estimated based on first principles:

$f'(x)=\frac{f(x+h)-f(x)}{h}$

Derivatives of continuous functions¶

At a desired $x$ for $f'(x)$ :

Choose an arbitrary $h$
Calculate derivative via first principles
Shrink $h$ and recalculate derivative
If the answer is drastically different, repeat step 3

Derivatives of discrete functions¶

Guesses are made based on the average slope between two points.

$f'(x_i)=\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}$

Divided differences¶

Using the next term, or a $\Delta x > 0$ indicates a forward divided difference (FDD).
Using the previous term, or a $\Delta x < 0$ indicates a backward divided difference (BDD).

The central divided difference averages both if $h$ or $\Delta x$ of the forward and backward DDs are equal.

$f'(x)=CDD=\frac{f(x+h)-f(x-h)}{2h}$

Higher order derivatives¶

Taking the Taylor expansion of the function or discrete set and then expanding it as necessary can return any order of derivative. This also applies for $x-h$ if positive and negative are alternated.

$f(x+h)=f(x)+f'(x)h+\frac{f''(x)}{2!}h^2+\frac{f'''(x)}{3!}h^3$

Example

To find second order derivatives:

$\begin{align*} f''(x)&=\frac{2f(x+h)-2f(x)-2f'(x)h}{h^2} \\ &=\frac{2f(x+h)-2f(x)-(f(x+h)-f(x-h))}{h^2} \\ &=\frac{f(x+h)-2f(x)+f(x-h)}{h^2} \end{align*}$

Example

$f''(3)$ if $f(x)=2e^{1.5x}$ and $h=0.1$ :

$\begin{align*} f''(3)&=\frac{f(3.1)-2\times2f(3)+f(2.9)}{0.1^2} \\ &=405.08 \end{align*}$

For discrete data:

If the desired point does not exist, differentiating the surrounding points to create a polynomial interpolation of the derivative may be close enough.

Example

t	0	10	15	20	22.5	30
v(t)	0	227.04	362.78	517.35	602.47	901.67

$v'(16)$ with FDD:

Using points $t=15,t=20$ :

$\begin{align*} v'(x)&=\frac{f(x+h)-f(x)}{h} \\ &=\frac{f(15+5)-f(15)}{5} \\ &=\frac{517.35-362.78}{5} \\ &=30.914 \end{align*}$

$v'(16)$ with Newton's first-order interpolation:

$\begin{align*} v(t)&=v(15)+\frac{v(20)-v(15)}{20-15}(t-15) \\ &=362.78+30.914(t-15) \\ &=-100.93+30.914t \\ v'(t)&=\frac{v(t+h)-v(t)}{2h} \\ &=\frac{v(16.1)-v(15.9)}{0.2} \\ &=30.914 \end{align*}$

If the spacing is not equal (to make DD impossible), again creating an interpolation may be close enough.
If data is noisy, regressing and then solving reduces random error.

Integrals¶

If you can represent a function as an $n$ -th order polynomial, you can approximate the integral with the integral of that polynomial.

Trapezoidal rule¶

The trapezoidal rule looks at the first order polynomial and

From $a$ to $b$ , if there are $n$ trapezoidal segments, where $h=\frac{b-a}{n}$ is the width of each segment:

$\int^b_af(x)dx=\frac{b-a}{2n}[f(a)+2(\sum^{n-1}_{i=1}f(a+ih))+f(b)]$

The error for the $i$ th trapezoidal segment is $|E_i|=\left|\frac{h^3}{12}\right|f''(x)$ . This can be approximated with a maximum value of $f''$ :

$\boxed{|E_T|\leq(b-a)\frac{h^2}{12}M}$

Simpson's 1/3 rule¶

This uses the second-order polynomial with two segments. Three points are usually used: $a,\frac{a+b}{2},b$ . Thus for two segments:

$\int^b_af(x)dx\approx\frac h 3\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$

For an arbitrary number of segments, as long as there are an even number of equal segments: $\int^b_af(x)dx=\frac{b-a}{3n}\left[f(x_0)+4\sum^{n-1}_{\substack{i=1 \\ \text{i is odd}}}f(x_i)+2\sum^{n-2}_{\substack{i=2 \\ \text{i is even}}}f(x_i)+f(x_n)\right]$

The error is: $|E_T|=(b-a)\frac{h^4}{180}M$

Ordinary differential equations¶

Initial value problems¶

These problems only have results for one value of $x$ .

Euler's method converts the function to the form $f(x,y)$ , where $y'=f(x,y)$ .

Example

$y'+2y=1.3e^{-x},y(0)=5\implies f(x,y)=1.3e^{-x}-2y,y(0)=5$

Where $h$ is the width of each estimation (lower is better):

$y_{n+1}=y_n+hf(x_n,y_n)$

Example

If $f(x,y)=2xy,h=0.1$ , $y_{n+1}=y_n+h2x_ny_n$

$y(1.1)=y(1)+0.1×2×1×\underbrace{y(1)}_{1 via IVP}=1.2 \\ y(1.2)=y(1.1)+0.1×2×1.1×\underbrace{y(1.1)}_{1.2}=1.464$

Heun's method uses Euler's formula as a predictor. Where $y^*$ is the Euler solution:

$y_{n+1}=y_n+h\frac{f(x_n,y_n)+f(x_{n+1},y^*_{n+1}}{2}$

Example

For $f(x,y)=2xy,h=0.1, y(1)=1$ :

Euler's formula returns $y^*_{n+1}=y_n+2hx_ny_n\implies y^*(1.1)=1.2$ .

Applying Heun's correction:

$\begin{align*} y(1.1)&=y(1)=0.1\frac{2×1×y(1)+2×1.1×y^*(1.1)}{2} \\ &=1+0.1\frac{2×1×1+2×1.1×1.2}{2} \\ &=1.232 \end{align*}$

The Runge-Kutta fourth-order method is the most accurate of the three methods:

$y_n+1=y_n+\frac 1 6(k_1+2k_2+2k_3+k_4)$

$k_1=hf(x_n,y_n)$
$k_2=hf(x_n+\tfrac 1 2h,y_n+\tfrac 1 2k_1)$
$k_3=hf(x_n+\tfrac 1 2 h, y_n+\tfrac 1 2 k_2)$
$k_4=hf(x_n+h,y_n+k_3)$

Higher order ODEs¶

Higher order ODEs can be solved by reducing them to first order ODEs by creating a system of equations. For a second order ODE: Let $y'=u$ .

$y'=u \\ u'=f(x,y,u)$

For each ODE, the any method can be used:

$y_{n+1}=y_n+hu_n \\ u_{n+1}=u_n+hf(x_n,y_n,u_n)$

Example

For $y''+xy'+y=0,y(0)=1,y'(0)=2,h=0.1$ :

$\begin{align*} y'&= u \\ u'&=-xu-y \\ y_1&=y_0+0.1u_0 \\ &=1+0.1×2 \\ &=1.2 \\ \\ u_1&=u_0+0.1×f(x_0,y_0,u_0) \\ &=u_0+0.1(-x_0u_0-y_0] \\ &=2+0.1(-0×2-1) \\ &=1.9 \end{align*}$

Boundary value problems¶

The finite difference method divides the interval between the boundary into $n$ sub-intervals, replacing derivatives with their first principles representations. Solving each $n-1$ equation returns a proper system of equations.

Example

For $y''+2y'+y=x^2, y(0)=0.2,y(1)=0.8,n=4\implies h=0.25$ : $x_0=0,x_1=0.25,x_2=0.5,x_3=0.75,x_4=1$

Replace with first principles:

$\frac{y_{i+1}-2y_i+y_{i-1}{h^2}+2\frac{y_{i+1}-y_i}{h}+y_i=x_i^2$