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

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.

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

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

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.