Functions

horner(c, x)

Evaluate a polynomial whose coefficients are given in ascending order in c, at the point x, using Horner's rule.

backsub(U,b)

Solve the upper-triangular linear system with matrix U and right-hand side vector b.

forwardsub(L,b)

Solve the lower-triangular linear system with matrix L and right-hand side vector b.

lufact(A)

Compute the LU factorization of square matrix A, returning the factors.

plufact(A)

Compute the PLU factorization of square matrix A, returning the triangular factors and a row permutation vector.

lsnormal(A, b)

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

lsqrfact(A, b)

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

qrfact(A)

QR factorization by Householder reflections. Returns Q and R.

fdjac(f, x₀ [,y₀])

Compute a finite-difference approximation of the Jacobian matrix for f at x₀, where y₀=f(x₀) may be given.

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.

newton(f, df_dx, x₁ [;maxiter,f tol, xtol])

Use Newton's method to find a root of f starting from x₁, where df_dx 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.

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.

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.

The optional keyword parameters set the maximum number of iterations and the stopping tolerance for values of f and changes in x.

fdweights(t, m)

Compute weights for the mth derivative of a function at zero using values at the nodes in vector t.

hatfun(t, k)

Create a piecewise linear hat function, where t is a vector of n+1 interpolation nodes and k is an integer in 0:n giving the index of the node where the hat function equals one.

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.

plinterp(t, y)

Construct a piecewise linear interpolating function for data values in y given at nodes in t.

spinterp(t, y)

Construct a cubic not-a-knot spline interpolating function for data values in y given at nodes in t.

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.

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.

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.

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.

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.

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.

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.

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.

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.)

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.

poisson(n)

Construct the finite-difference Laplacian matrix for an n by n interior lattice.

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.

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.

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.

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.

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.

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.

polyinterp(t, y)

Construct a callable polynomial interpolant through the points in vectors t, y using the barycentric interpolation formula.

triginterp(t, y)

Construct the trigonometric interpolant for the points defined by vectors t and y.

bvp(ϕ, xspan, g₁, g₂, 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.

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.

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.

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.

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.

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'.

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.

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.

Evaluate Chebyshev interpolant with nodes x, values v, at point ξ

elliptic(ϕ, g, m, xspan, n, yspan)

Solve the elliptic PDE ϕ(x, y, u, ux, uxx, uy, uyy) = 0 on the rectangle xspanxyspan, 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.

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.

tensorgrid(x, y)

Create a tensor grid for a rectangle from its 1d projections x and y. Returns unvec to reshape a vector into a 2d array, mtx to evaluate a function on the grid, X, Y to give the grid coordinates, and boolean array is_boundary to identify the boundary points.