Scientific computing develops numerical methods for problems whose exact solutions are unavailable, inconvenient, or too expensive to evaluate. The central workflow is:

1. formulate a mathematical model;
2. replace derivatives, integrals, or functions by computable approximations;
3. derive an error estimate;
4. study convergence and stability; and
5. validate the method with computation.

This note collects the interpolation, numerical integration, initial-value problem, Runge–Kutta, stability, and nonlinear-pendulum material from the course.

1 Polynomial Interpolation link

Suppose that a function is known at distinct nodes

We seek a polynomial $P_n$ of degree at most $n$ satisfying

1.1 Lagrange Basis link

For each node $x_i$, construct a polynomial that vanishes at every other node. Its unnormalized numerator is

At $x=x_i$, this numerator has the value

Dividing by this value normalizes the polynomial to one:

Indeed, if $k\ne i$, one factor is $x_k-x_k=0$, so $L_i(x_k)=0$. If $k=i$, every numerator and denominator factor is equal, so $L_i(x_i)=1$. Therefore,

Taking a linear combination with weights $y_i$ gives the Lagrange interpolating polynomial:

At a node $x_k$,

so the polynomial interpolates every data point.

1.2 Example link

For $(x_0,y_0)=(1,2)$, $(x_1,y_1)=(2,4)$, and $(x_2,y_2)=(4,8)$, the basis function associated with $x_1=2$ is

At $x=3$,

Consequently,

1.3 Uniqueness link

Suppose $P_n$ and $Q_n$ are two polynomials of degree at most $n$ that interpolate the same $n+1$ data points. Let

Then $R(x_i)=0$ for every $i$, so $R$ has at least $n+1$ distinct roots. A nonzero polynomial of degree at most $n$ can have at most $n$ roots. Hence $R$ is the zero polynomial and

The interpolating polynomial is therefore unique.

1.4 Partition of Unity and Affine Invariance link

Apply interpolation to the constant function $f(x)=1$. Since $f(x_i)=1$,

Thus the basis functions form a partition of unity: their weights always sum to one. In particular, interpolation preserves constants and affine combinations.

If the independent variable is changed by an affine map $x=at+b$, the same construction in the $t$-coordinate produces the correspondingly transformed interpolating polynomial. Translation and scaling change the coordinates, but not the interpolation structure.

1.5 Nodal Polynomial Form link

Define the nodal polynomial

Removing the factor corresponding to $x_i$ gives

Differentiating the product and evaluating at $x_i$, every term vanishes except the term in which the factor $(x-x_i)$ is differentiated:

Therefore,

and the interpolant can also be written as

1.6 Interpolation Error link

Assume $f\in C^{n+1}[a,b]$, with all nodes and the evaluation point $x$ in the interval. Define

At the interpolation nodes, $E(x_i)=0$. For a point $x$ not equal to a node, choose a coefficient $c(x)$ so that

Now define an auxiliary function of $z$:

By construction, $F(x_i)=0$ for $i=0,\ldots,n$, and $F(x)=0$. Thus $F$ has $n+2$ distinct roots. Repeated use of Rolle’s theorem gives a point $\xi$ between the nodes and $x$ such that

Since $P_n^{(n+1)}=0$ and $Q^{(n+1)}=(n+1)!$,

Therefore,

and the error formula is

If $f$ is a polynomial of degree $k\le n$, then $f^{(n+1)}=0$. Hence the interpolation is exact: $E(x)=0$ for every $x$.

1.7 Derivatives of the Interpolant link

Differentiate the interpolant term by term:

Starting with

differentiate once:

At another node $x_j$, $j\ne i$, $L_i(x_j)=0$, so

At $x=x_i$, differentiate the identity once more:

Evaluating at $x_i$ gives

The derivative error follows by differentiating the interpolation error. Writing the intermediate point as $\xi=\xi(x)$,

At a node $x_i$, the second term disappears because $Q(x_i)=0$. Hence

1.8 Runge’s Phenomenon link

Increasing the number of nodes does not always improve interpolation. For

use equally spaced nodes $x_i=-1+2i/n$. The unique polynomial is

For large $n$, the product factors become badly conditioned near the endpoints and the interpolant can oscillate strongly there. This is Runge’s phenomenon. More nodes are not sufficient by themselves; node placement, such as Chebyshev nodes, is also important.

2 Numerical Integration link

Let

and divide $[a,b]$ into $n$ equal subintervals. Set

2.1 Left and Right Rectangle Rules link

The left-point approximation is

and the right-point approximation is

For the right rule, write the error as a sum of local errors:

Taylor’s theorem about the right endpoint gives

Therefore,

If $M_1=\max_{[a,b]}|f'(x)|$, the leading term satisfies

Summing $n$ intervals,

The same argument about the left endpoint gives

Thus both endpoint rules are first-order methods: their global error is $O(h)=O(n^{-1})$.

2.2 Midpoint Rule link

Let $m_i=(x_{i-1}+x_i)/2$. Approximate $f$ about the midpoint:

The local midpoint error is

The linear term cancels by symmetry:

For the quadratic term, use $t=x-m_i$, with limits $-h/2$ and $h/2$:

If $M_2=\max|f''|$, then

The midpoint rule is second-order.

2.3 Trapezoidal Rule link

On $[x_{i-1},x_i]$, the trapezoidal approximation is

Summing gives

The local trapezoidal error is

so

The trapezoidal rule is also second-order.

2.4 Simpson’s Rule from Lagrange Interpolation link

Use three consecutive points $x_{i-1}$, $x_i$, and $x_{i+1}$, where the spacing is $h$. The quadratic Lagrange interpolant is

where

Set $x=x_i+ht$, so $dx=h\,dt$ and $t\in[-1,1]$. Then

because

Similarly,

because

and

Therefore,

Adding the $n/2$ panels, with $n$ even, gives the composite Simpson rule:

Taylor expansion about the midpoint shows that the odd powers cancel. The first remaining neglected term is proportional to $h^5f^{(4)}$, so the global error is $O(h^4)$.

2.5 Numerical Experiment link

For a smooth test function such as $f(x)=e^{-x^2}$, compute a high-precision reference value and compare the five rules for $N=10,20,40,80,160$. If $E_N$ denotes the error, the observed order is estimated by

The expected orders are one for left/right rules, two for midpoint and trapezoidal rules, and four for Simpson’s rule. Symmetry can occasionally produce a higher observed order for a particular function, but the general theoretical guarantees remain those listed above.

3 Initial-Value Problems link

Consider

Integrating from $t$ to $t+h$ and applying the fundamental theorem of calculus gives the exact identity

Every one-step method below approximates this integral in a different way.

3.1 Euler Method link

Approximate the integrand by its left endpoint:

Then

and, with $t_n=nh$ and $v_n\approx u(t_n)$,

Writing $K_1=f(t_n,v_n)$ gives $v_{n+1}=v_n+hK_1$.

3.2 Midpoint Method link

Approximate the integral with the slope at the midpoint:

The midpoint value is unknown, so predict it by a half Euler step:

Thus

3.3 Heun’s Method link

Apply the trapezoidal rule to the IVP integral:

Predict the unknown endpoint by Euler:

Then

Both midpoint and Heun methods have local error $O(h^3)$ and global error $O(h^2)$.

3.4 Kutta’s Third-Order Method link

One third-order method from the course is

The weights are chosen so that the Taylor expansion agrees through terms of order three.

3.5 Classical RK4 link

Use two midpoint estimates, one based on $K_1$ and one based on $K_2$:

The update is

The four slopes approximate the integral with a Simpson-like weighted average. Taylor matching gives local truncation error $O(h^5)$ and global error $O(h^4)$.

4 Error, Consistency, and Stability link

4.1 Euler Residual link

Let $u_n=u(t_n)$ be the exact solution and $v_n$ the numerical solution. The residual of Euler’s scheme is

Taylor’s theorem gives

so, because $u'=f(t,u)$,

Since $r_n\to0$ as $h\to0$, the method is consistent and first-order accurate as an approximation to the differential equation.

4.2 Euler Error Equation link

Let $e_n=u_n-v_n$. Subtract the numerical update from the exact equation with residual:

Then

By the mean-value theorem,

where $\eta_n$ lies between $u_n$ and $v_n$. Hence

Assume $|u''(t)|\le M$ and $|f_u(t,u)|\le L$. Then

Iterating from $e_0=0$,

Since $1+hL\le e^{hL}$ and $T=(n+1)h$,

For fixed $T$, this tends to zero linearly with $h$.

4.3 Stability and Amplification Factors link

For the test equation

Euler gives

The amplification factor is

Absolute stability requires

For $\lambda<0$, this becomes

so

For example, $u'=-5u$ requires $h\le0.4$ for stable explicit Euler integration. For a general Runge–Kutta method, the stability condition is $|R(h\lambda)|\le1$; for a system, $h\lambda_i$ must lie in the method’s stability region for every relevant Jacobian eigenvalue $\lambda_i$.

5 Nonlinear Simple Pendulum link

5.1 Physical Model link

Consider a point mass $m$ attached to a massless, inextensible string of length $L$ in a uniform gravitational field $g$. Let $\theta(t)$ be the angle measured from the downward vertical. The arc length is

The gravitational force has tangential component $-mg\sin\theta$. Newton’s second law in the tangential direction gives

Since $s=L\theta$,

Substituting and cancelling $m$ yields

The initial-value problem is

5.2 Small-Angle Solution link

For small angles measured in radians,

The nonlinear equation becomes

Let $\omega_* =\sqrt{g/L}$. The general solution is

At $t=0$,

Differentiating,

so

Therefore,

Writing this as $C\cos(\omega_*t-\phi)$ gives

The period satisfies $\omega_*T=2\pi$, hence

This small-angle period is independent of amplitude and initial angular velocity.

5.3 Large-Angle Period link

For large angles, retain $\sin\theta$. Conservation of energy gives, for release from rest at amplitude $\theta_0$,

Using $1-\cos\theta=2\sin^2(\theta/2)$,

Set $k=\sin(\theta_0/2)$ and use the substitution

After simplifying the quarter-period integral,

Thus

where $K(k)$ is the complete elliptic integral of the first kind. In software such as SciPy, ellipk commonly takes the parameter $m=k^2$:

  from scipy.special import ellipk
import numpy as np

def exact_period(L, g, theta0):
    m = np.sin(theta0 / 2.0) ** 2
    return 4.0 * np.sqrt(L / g) * ellipk(m)
  

5.4 First-Order System link

Define angular velocity $\omega=\theta'$. Then

RK4 is applied to this vector equation component by component.

5.5 RK4 Pendulum Simulation link

  import numpy as np

def pendulum_derivatives(t, y, g, L):
    theta, omega = y
    return np.array([omega, -(g / L) * np.sin(theta)])

def rk4_step(f, t, y, h, g, L):
    k1 = f(t, y, g, L)
    k2 = f(t + h / 2, y + h * k1 / 2, g, L)
    k3 = f(t + h / 2, y + h * k2 / 2, g, L)
    k4 = f(t + h, y + h * k3, g, L)
    return y + h * (k1 + 2*k2 + 2*k3 + k4) / 6

def simulate(theta0=np.pi / 4, L=1.0, g=9.81,
             t_end=10.0, h=0.05):
    times = np.arange(0.0, t_end, h)
    states = np.zeros((len(times), 2))
    states[0] = [theta0, 0.0]
    for i in range(1, len(times)):
        states[i] = rk4_step(
            pendulum_derivatives, times[i - 1], states[i - 1], h, g, L
        )
    return times, states
  

The trajectory can be plotted as $\theta(t)$, and the phase portrait is obtained by plotting angular velocity against angle.

5.6 Exact Trajectory and Convergence link

For release from rest, the nonlinear trajectory can be expressed with Jacobi elliptic functions. With

one representation is

The numerical error at the time nodes is

Repeat the experiment with several step sizes, for example

If RK4 is in its asymptotic regime, halving $h$ reduces the global error by approximately $2^4=16$. A log-log graph of error versus $h$ has slope approximately four:

5.7 Period Extraction by Zero Crossings link

A numerical period can be obtained from successive turning points, where $\omega(t)=0$. If $\omega_n$ and $\omega_{n+1}$ have opposite signs, linearly interpolate the zero:

The time between every second turning point is one full oscillation period. Comparing this value with the elliptic-integral period measures the combined integration and event-location error.

5.8 Symplectic Euler link

The pendulum is conservative, so long-time energy behavior matters. Explicit RK4 is accurate over a fixed interval but may slowly distort energy over very long integrations. A semi-implicit Euler method updates velocity first and then angle:

It is only first-order accurate, but its geometric structure often prevents the systematic energy drift seen in ordinary explicit Euler.

6 Higher-Order Runge–Kutta Methods link

An explicit Runge–Kutta method obtains several stage slopes and combines them with weights. Increasing the order requires satisfying increasingly many nonlinear order conditions. These conditions can be organized using rooted trees and B-series.

For explicit methods of order $p\ge5$, the Butcher barrier implies that more than $p$ stages are required. For example, fifth order needs at least six stages and sixth order needs at least seven stages. The additional evaluations reduce the efficiency advantage of a fixed high-order method.

Embedded methods such as RKF45 and Dormand–Prince use the same stages with two weight vectors $b_i$ and $b_i^*$:

The difference

estimates the local error. The solver can then increase $h$ when the estimate is small and decrease $h$ when it is large. This produces adaptive step-size control without evaluating the differential equation twice with two unrelated methods.

7 Practical Checklist link

• Check that interpolation nodes are distinct.
• Use the interpolation error formula to identify smoothness and node-placement effects.
• Distinguish local truncation error from accumulated global error.
• Verify an observed convergence order with a refinement table or log-log plot.
• For ODEs, check both accuracy and stability; a small local error does not guarantee a useful long-time solution.
• For conservative dynamics, monitor energy and phase-space behavior.
• Keep angle units in radians when evaluating trigonometric functions.
• When using scipy.special.ellipk, pass the parameter $m=k^2$, not the modulus $k$ itself.