# Functions

## Chapter 1: Introduction

`FundamentalsNumericalComputation.horner`

— Method`horner(c,x)`

Evaluate a polynomial whose coefficients are given in ascending order in `c`

, at the point `x`

, using Horner's rule.

## Chapter 2: Square linear systems

`FundamentalsNumericalComputation.backsub`

— Method`backsub(U,b)`

Solve the upper-triangular linear system with matrix `U`

and right-hand side vector `b`

.

`FundamentalsNumericalComputation.forwardsub`

— Method`forwardsub(L,b)`

Solve the lower-triangular linear system with matrix `L`

and right-hand side vector `b`

.

`FundamentalsNumericalComputation.lufact`

— Method`lufact(A)`

Compute the LU factorization of square matrix `A`

, returning the factors.

`FundamentalsNumericalComputation.plufact`

— Method`plufact(A)`

Compute the PLU factorization of square matrix `A`

, returning the triangular factors and a row permutation vector.

## Chapter 3: Overdetermined linear systems

`FundamentalsNumericalComputation.lsnormal`

— Method`lsnormal(A,b)`

Solve a linear least-squares problem by the normal equations. Returns the minimizer of ||b-Ax||.

`FundamentalsNumericalComputation.lsqrfact`

— Method`lsqrfact(A,b)`

Solve a linear least-squares problem by QR factorization. Returns the minimizer of ||b-Ax||.

`FundamentalsNumericalComputation.qrfact`

— Method`qrfact(A)`

QR factorization by Householder reflections. Returns Q and R.

## Chapter 4: Roots of nonlinear equations

`FundamentalsNumericalComputation.fdjac`

— Function`fdjac(f,x₀[,y₀])`

Compute a finite-difference approximation of the Jacobian matrix for `f`

at `x₀`

, where `y₀`

=`f(x₀)`

may be given.

`FundamentalsNumericalComputation.levenberg`

— Method`levenberg(f,x₁[;maxiter,ftol,xtol])`

Use Levenberg's quasi-Newton iteration to find a root of the system `f`

starting from `x₁`

Returns the history of root estimates as a vector of vectors.

The optional keyword parameters set the maximum number of iterations and the stopping tolerance for values of `f`

and changes in `x`

.

`FundamentalsNumericalComputation.newton`

— Method`newton(f,dfdx,x₁[;maxiter,ftol,xtol])`

Use Newton's method to find a root of `f`

starting from `x₁`

, where `dfdx`

is the derivative of `f`

. Returns a vector of root estimates.

The optional keyword parameters set the maximum number of iterations and the stopping tolerance for values of `f`

and changes in `x`

.

`FundamentalsNumericalComputation.newtonsys`

— Method`newtonsys(f,jac,x₁[;maxiter,ftol,xtol])`

Use Newton's method to find a root of a system of equations, starting from `x₁`

. The functions `f`

and `jac`

should return the residual vector and the Jacobian matrix, respectively. Returns the history of root estimates as a vector of vectors.

The optional keyword parameters set the maximum number of iterations and the stopping tolerance for values of `f`

and changes in `x`

.

`FundamentalsNumericalComputation.secant`

— Method`secant(f,x₁,x₂[;maxiter,ftol,xtol])`

Use the secant method to find a root of `f`

starting from `x₁`

and `x₂`

. Returns a vector of root estimates.

`f`

and changes in `x`

.

## Chapter 5: Piecewise interpolation

`FundamentalsNumericalComputation.fdweights`

— Method`fdweights(t,m)`

Compute weights for the `m`

th derivative of a function at zero using values at the nodes in vector `t`

.

`FundamentalsNumericalComputation.intadapt`

— Function`intadapt(f,a,b,tol)`

Adaptively integrate `f`

over [`a`

,`b`

] to within target error tolerance `tol`

. Returns the estimate and a vector of evaluation nodes.

`FundamentalsNumericalComputation.plinterp`

— Method`plinterp(t,y)`

Construct a piecewise linear interpolating function for data values in `y`

given at nodes in `t`

.

`FundamentalsNumericalComputation.spinterp`

— Method`spinterp(t,y)`

Construct a cubic not-a-knot spline interpolating function for data values in `y`

given at nodes in `t`

.

`FundamentalsNumericalComputation.trapezoid`

— Method`trapezoid(f,a,b,n)`

Apply the trapezoid integration formula for integrand `f`

over interval [`a`

,`b`

], broken up into `n`

equal pieces. Returns the estimate, a vector of nodes, and a vector of integrand values at the nodes.

## Chapter 6: Initial-value problems for ODEs

`FundamentalsNumericalComputation.ab4`

— Method`ab4(ivp,n)`

Apply the Adams-Bashforth 4th order method to solve the given IVP using `n`

time steps. Returns a vector of times and a vector of solution values.

`FundamentalsNumericalComputation.am2`

— Method`am2(ivp,n)`

Apply the Adams-Moulton 2nd order method to solve given IVP using `n`

time steps. Returns a vector of times and a vector of solution values.

`FundamentalsNumericalComputation.euler`

— Method`euler(ivp,n)`

Apply Euler's method to solve the given IVP using `n`

time steps. Returns a vector of times and a vector of solution values.

`FundamentalsNumericalComputation.ie2`

— Method`ie2(ivp,n)`

Apply the Improved Euler method to solve the given IVP using `n`

time steps. Returns a vector of times and a vector of solution values.

`FundamentalsNumericalComputation.rk23`

— Method`rk23(ivp,tol)`

Apply an adaptive embedded RK formula pair to solve given IVP with estimated error `tol`

. Returns a vector of times and a vector of solution values.

`FundamentalsNumericalComputation.rk4`

— Method`rk4(ivp,n)`

Apply the common Runge-Kutta 4th order method to solve the given IVP using `n`

time steps. Returns a vector of times and a vector of solution values.

## Chapter 7: Matrix analysis

## Chapter 8: Krylov methods in linear algebra

`FundamentalsNumericalComputation.arnoldi`

— Method`arnoldi(A,u,m)`

Perform the Arnoldi iteration for `A`

starting with vector `u`

, out to the Krylov subspace of degree `m`

. Returns the orthonormal basis (`m`

+1 columns) and the upper Hessenberg `H`

of size `m`

+1 by `m`

.

`FundamentalsNumericalComputation.gmres`

— Method`gmres(A,b,m)`

Do `m`

iterations of GMRES for the linear system `A`

*x=`b`

. Returns the final solution estimate x and a vector with the history of residual norms. (This function is for demo only, not practical use.)

`FundamentalsNumericalComputation.inviter`

— Method`inviter(A,s,numiter)`

Perform `numiter`

inverse iterations with the matrix `A`

and shift `s`

, starting from a random vector. Returns a vector of eigenvalue estimates and the final eigenvector approximation.

`FundamentalsNumericalComputation.poisson`

— Method`poisson(n)`

Construct the finite-difference Laplacian matrix for an `n`

by `n`

interior lattice.

`FundamentalsNumericalComputation.poweriter`

— Method`poweriter(A,numiter)`

Perform `numiter`

power iterations with the matrix `A`

, starting from a random vector. Returns a vector of eigenvalue estimates and the final eigenvector approximation.

`FundamentalsNumericalComputation.sprandsym`

— Method```
sprandsym(n,density,λ)
sprandsym(n,density,rcond)
```

Construct a randomized `n`

by `n`

symmetric sparse matrix of approximate density `density`

. For vector `λ`

, the matrix has eigenvalues as prescribed by λ. For scalar `rcond`

, the matrix has condition number equal to 1/`rcond`

.

## Chapter 9: Global function approximation

`FundamentalsNumericalComputation.ccint`

— Method`ccint(f,n)`

Perform Clenshaw-Curtis integration for the function `f`

on `n`

+1 nodes in [-1,1]. Returns the integral estimate and a vector of the nodes used. Note: `n`

must be even.

`FundamentalsNumericalComputation.glint`

— Method`glint(f,n)`

Perform Gauss-Legendre integration for the function `f`

on `n`

nodes in (-1,1). Returns the integral estimate and a vector of the nodes used.

`FundamentalsNumericalComputation.intinf`

— Method`intinf(f,tol)`

Perform adaptive doubly-exponential integration of function `f`

over (-Inf,Inf), with error tolerance `tol`

. Returns the integral estimate and a vector of the nodes used.

`FundamentalsNumericalComputation.intsing`

— Method`intsing(f,tol)`

Adaptively integrate function `f`

over (0,1), where `f`

may be singular at zero, with error tolerance `tol`

. Returns the integral estimate and a vector of the nodes used.

`FundamentalsNumericalComputation.polyinterp`

— Method`polyinterp(t,y)`

Construct a callable polynomial interpolant through the points in vectors `t`

,`y`

using the barycentric interpolation formula.

`FundamentalsNumericalComputation.triginterp`

— Method`triginterp(t,y)`

Construct the trigonometric interpolant for the points defined by vectors `t`

and `y`

.

## Chapter 10: Boundary-value problems

`FundamentalsNumericalComputation.bvp`

— Method`bvp(ϕ,xspan,lval,lder,rval,rder,init)`

Finite differences to solve a two-point boundary value problem with ODE u'' = `ϕ`

(x,u,u') for x in `xspan`

, left boundary condition `g₁`

(u,u')=0, and right boundary condition `g₂`

(u,u')=0. The value `init`

is an initial estimate for the values of the solution u at equally spaced values of x, which also sets the number of nodes.

Returns vectors for the nodes and the values of u.

`FundamentalsNumericalComputation.bvplin`

— Method`bvplin(p,q,r,xspan,lval,rval,n)`

Use finite differences to solve a linear bopundary value problem. The ODE is u''+`p`

(x)u'+`q`

(x)u = `r`

(x) on the interval `xspan`

, with endpoint function values given as `lval`

and `rval`

. There will be `n`

+1 equally spaced nodes, including the endpoints.

Returns vectors of the nodes and the solution values.

`FundamentalsNumericalComputation.diffcheb`

— Method`diffcheb(n,xspan)`

Compute Chebyshev differentiation matrices on `n`

+1 points in the interval `xspan`

. Returns a vector of nodes and the matrices for the first and second derivatives.

`FundamentalsNumericalComputation.diffmat2`

— Method`diffmat2(n,xspan)`

Compute 2nd-order-accurate differentiation matrices on `n`

+1 points in the interval `xspan`

. Returns a vector of nodes and the matrices for the first and second derivatives.

`FundamentalsNumericalComputation.fem`

— Method`fem(c,s,f,a,b,n)`

Use a piecewise linear finite element method to solve a two-point boundary value problem. The ODE is (`c`

(x)u')' + `s`

(x)u = `f`

(x) on the interval [`a`

,`b`

], and the boundary values are zero. The discretization uses `n`

equal subintervals.

Return vectors for the nodes and the values of u.

`FundamentalsNumericalComputation.shoot`

— Function`shoot(ϕ,xspan,g₁,g₂,init)`

Shooting method to solve a two-point boundary value problem with ODE u'' = `ϕ`

(x,u,u') for x in `xspan`

, left boundary condition `g₁`

(u,u')=0, and right boundary condition `g₂`

(u,u')=0. The value `init`

is an initial estimate for vector [u,u'] at x=a.

Returns vectors for the nodes, the solution u, and derivative u'.

## Chapter 11: Diffusion equations

`FundamentalsNumericalComputation.diffper`

— Method`diffper(n,xspan)`

Construct 2nd-order differentiation matrices for functions with periodic end conditions, using `n`

unique nodes in the interval `xspan`

. Returns a vector of nodes and the matrices for the first and second derivatives.

`FundamentalsNumericalComputation.parabolic`

— Method`parabolic(ϕ,xspan,m,g₁,g₂,tspan,init)`

Solve a parabolic PDE by the method of lines. The PDE is ∂u/∂t = `ϕ`

(t,x,u,∂u/∂x,∂^2u/∂x^2), `xspan`

gives the space domain, m gives the degree of a Chebyshev spectral discretization, `g₁`

and `g₂`

are functions of (u,∂u/∂x) at the domain ends that should be made zero, `tspan`

is the time domain, and `init`

is a function of x that gives the initial condition. Returns a vector `x`

and a function of t that gives the semidiscrete solution at `x`

.

## Chapter 12: Advection equations

## Chapter 13: Two-dimensional problems

`FundamentalsNumericalComputation.chebinterp`

— MethodEvaluate Chebyshev interpolant with nodes x, values v, at point ξ

`FundamentalsNumericalComputation.elliptic`

— Method`elliptic(ϕ,g,m,xspan,n,yspan)`

Solve the elliptic PDE `ϕ`

(x,y,u,u*x,u*xx,u*y,u*yy)=0 on the rectangle `xspan`

x `yspan`

, subject to `g`

(x,y)=0 on the boundary. Uses `m`

+1 points in x by `n`

+1 points in y in a Chebyshev discretization.

Returns vectors defining the grid and a matrix of grid solution values.

`FundamentalsNumericalComputation.poissonfd`

— Method`poissonfd(f,g,m,xspan,n,yspan)`

Solve Poisson's equation on a rectangle by finite differences. Function `f`

is the forcing function and function `g`

gives the Dirichlet boundary condition. The rectangle is the tensor product of intervals `xspan`

and `yspan`

, and the discretization uses `m`

+1 and `n`

+1 points in the two coordinates.

Returns vectors defining the grid and a matrix of grid solution values.