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.

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,x0,y0)

Compute a finite-difference approximation of the Jacobian matrix for f at x0, where y0=f(x0) is given.

levenberg(f,x₁,tol)

Use Levenberg's quasi-Newton iteration to find a root of the system f, starting from x₁, with tol as the stopping tolerance in both step size and residual norm. Returns the history of root estimates as a vector of vectors.

newton(f,dfdx,x₁)

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.

newtonsys(f,jac,x₁)

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.

secant(f,x₁,x₂)

Use the secant method to find a root of f starting from x₁ and x₂. Returns a vector of root estimates.

fdweights(t,m)

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

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.

intde(f,h,M)

Perform doubly-exponential integration of function f over (-Inf,Inf), using discretization size h and truncation point M. Returns the integral estimate and a vector of the nodes used.

intsing(f,h,δ)

Integrate function f (possibly singular at 1 and -1) over [-1+δ,1-δ] using discretization size h. 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(phi,xspan,lval,lder,rval,rder,init)

Use finite differences to solve a two-point boundary value problem. The ODE is u'' = phi(x,u,u') for x in xspan. Specify a function value or derivative at the left endpoint using lval and lder, respectively, and similarly for the right endpoint using rval and rder. (Use an empty array to denote an unknown quantity.) The value init is an initial guess for whichever value is missing at the left endpoint.

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(phi,xspan,lval,lder,rval,rder,init)

Use the shooting method to solve a two-point boundary value problem. The ODE is u'' = phi(x,u,u') for x in xspan. Specify a function value or derivative at the left endpoint using lval and lder, respectively, and similarly for the right endpoint using rval and rder. (Use an empty array to denote an unknown quantity.) The value init is an initial guess for whichever value is missing at the left endpoint.

Returns vectors for the nodes, the values of u, and the values of 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.

ndgrid(x,y,...)

Given $d$ vector inputs, returns $d$` matrices representing the coordinate functions on the tensor product grid.

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 matrices of the grid solution values and the coordinate functions.

rectdisc(m,xspan,n,yspan)

Construct matrices and helpers for finite-difference discretization of a rectangle that is the tensor product of intervals xspan and yspan, using m+1 and n+1 points in the two coordinates.