GNU Octave is also freely redistributable software. You may redistribute
it and/or modify it under the terms of the GNU General Public License as
published by the
Free Software Foundation. The GPL is included in this manual in section
GNU GENERAL PUBLIC LICENSE.
This document corresponds to Octave version 2.1.33 and the modifications
made to it for Wayne State College, Nebraska.
octave:3> function y = funn(x)
> y = x.^3 +x.^2 -3*x -3;
> endfunction
octave:4>
A script-file can be written using any text editor such as vi, emacs, or gedit. The file must contain the same command statements as used in Octave's command window. The file is then saved in the present working directory using the name of the function with a ".m" suffix, i.e. "funn.m".
You can use the following shell commands in Octave.
Additional routines for use in WSC's Numerical Analysis course and similar courses, have been added to Octave. These routines are separated into four different classes.roots
bisect (bisection method for algebraic equations) secant (secant method for algebraic equations) newton (Newton's method for algebraic equations) muller (Muller's method for algebraic equations) bairstow (Bairstow's method for polynomial equations) rref (reduced row echelon form) jacobi (Jacobi iterative method) seidel (Gauss-Seidel iterative method) forward (forward substitution) back (backward substitution) tridia (Gaussian elimination on a tridiagonal matrix)
Under Construction
usage: [x, y] = bisect(FunFcn, a, b, tolerance [optional], maximum [optional])
bisect approximates the zeros to a function based on a sign change
between
f(a) and f(b). The function is described by
a m-file. The routine is
performed until either the maximum number of iterations has been
exhausted
or until an approximate solution has been found within the desired
tolerance,
i.e. |x(i-1) - x(i)| <= tolerance.
A table displaying the values for each
iteration is displayed on the screen. The x values and the
y
values are also
returned.
INPUT:
FunFcn - a string name of an m-file containing
an algebraic equation of
the form y = f(x). x is the independent variable and y
is the
returned dependent variable.
a - first bracketing
value for x.
b - second bracketing
value for x. a and b must be selected so that
f(a)*f(b)
< 0.
tolerance - bisect will perform iterations until the difference,
|x(i-1) - x(i)| <= tolerance. The default tolerance is
0.001.
[optional]
maximum - bisect will perform up to a maximum
number of iterations.
The default maximum is 100. [optional]
OUTPUT:
x - returned x
values.
y - returned y
values.
A table displaying each iteration is displayed on the screen.
usage: [x, y] = secant(FunFcn, a, b, tolerance [optional], maximum [optional])
secant approximates the zeros to a function. Unlike the bisect method,
it is
not based on a sign change between f(a) and f(b).
The function is described
by a m-file. The routine is performed until either the maximum
number of
iterations has been exhausted or until an approximate solution
has been found
within the desired tolerance, i.e. |x(i-1) - x(i)|
<= tolerance. A table
displaying the values for each iteration is displayed on the screen.
The x values and the y values are also returned.
INPUT:
FunFcn - a string name of an m-file containing
an algebraic equation of
the form y = f(x). x is the independent variable and y
is the
returned dependent variable.
a - first bracketing
value for x.
b - second bracketing
value for x. a and b should be selected so
that f(a)*f(b) < 0.
tolerance - secant will perform iterations until the difference,
|x(i-1) - x(i)| <= tolerance. The default tolerance is
0.001.
[optional]
maximum - secant will perform up to a maximum
number of iterations.
The default maximum is 100. [optional]
OUTPUT:
x - returned x
values.
y - returned y
values.
A table displaying each iteration is displayed on the screen.
usage: [x, y] = newton(FunFcn, FunFcn_prime, x1, tolerance [optional],
maximum [optional])
newton approximates the zeros to a function based on an initial
guess, x1,
to the solution. The function is described by a m-file and so too
is its
derivative. The routine is performed until either the maximum number
of
iterations has been exhausted or until an approximate solution
has been
found within the desired tolerance, i.e. |x(i-1) - x(i)|
<= tolerance.
A table displaying the values for each iteration is displayed on
the
screen. The x values and the y values are also returned.
INPUT:
FunFcn - a string name of an m-file containing
an algebraic equation of
the form y = f(x). x is the independent variable and y
is the
returned dependent variable.
FunFcn_prime - a string name of an m-file containing an algebraic
equation of
the form yprime = f(x). x is the independent variable and
y
is the
returned derivative of the original function.
x1 - the initial guess
at the solution.
tolerance - newton will perform iterations until the difference,
|x(i-1) - x(i)| <= tolerance. The default tolerance is
0.001.
[optional]
maximum - newton will perform up to a maximum
number of iterations.
The default maximum is 100. [optional]
OUTPUT:
x - returned x
values.
y - returned y
values.
A table displaying each iteration is displayed on the screen.
usage: [x, y] = muller(FunFcn, a, b, tolerance [optional], maximum [optional])
muller approximates the zeros to a function based on a sign change
between
f(a) and f(b). The function is described by
a m-file. The routine is
performed until either the maximum number of iterations has been
exhausted
or until an approximate solution has been found within the desired
tolerance,
i.e. |x(i-1) - x(i)| <= tolerance.
A table displaying the values for each
iteration is displayed on the screen. The x values and the
y
values are also
returned.
INPUT:
FunFcn - a string name of an m-file containing
an algebraic equation of
the form y = f(x). x is the independent variable and y
is the
returned dependent variable.
a - first bracketing
value for x.
b - second bracketing
value for x. a and b must be selected so that
f(a)*f(b)
< 0.
tolerance - muller will perform iterations until the difference,
|x(i-1) - x(i)| <= tolerance. The default tolerance is
0.001.
[optional]
maximum - muller will perform up to a maximum
number of iterations.
The default maximum is 100. [optional]
OUTPUT:
x - returned x
values.
y - returned y
values.
A table displaying each iteration is displayed on the screen.
usage: [rts, quot] = bairstow(polynomial, r_0 [optional], s_0 [optional],
maximum [optional], tolerance [optional])
bairstow approximates the zeros to a polynomial of degree n (n >
2).
The polynomial must be defined as a row vector of real coefficients.
The
polynomial is of the form x^n + a(1)x^(n-1) + a(2)x^(n-2) +
... + a(n).
The bairstow routine performs quadratic synthetic division. This
process
is performed until the approximate solution falls within the defined
tolerance or the maximum number of iterations have been performed.
A pair
of roots and the quotient are returned.
INPUT:
Polynomial - a row vector containing the real coefficients of a
polynomial,
i.e.[1, a(1), a(2), ... ,a(n)]. Gaps must be filled with zeros.
r_0 - the initial value for
the linear coefficient of the divisor.
Default value is 0.31. [optional]
s_0 - the initial value for
the constant of the divisor. Default value
is 0.71. [optional]
maximum - bairstow will perform up to a maximum
number of iterations.
The default maximum is 100. [optional]
tolerance - bairstow will perform iterations until the
(|dr| and |ds|) <= tolerance. The default tolerance is 0.001.
[optional]
OUTPUT:
rts - the approximate roots
for the quadratic divisor.
quot - the coeeficients to the (n-2)
degree quotient.
Reduced row echelon form.
usage: A = rref(A, tol [optional])
rref produces the reduced row echelon form of matrix A. The matrix A should be an augmented matrix of size m X m+1.
INPUT:
A - a matrix.
The matrix should be an augmented matrix of
dimensions m X m+1.
tol - rref will convert
matrix elements, A(i,j) < tol to 0. The
tolerance is set based on the infinite norm of the matrix.
OUTPUT:
A - the reduced
row echelon form of the matrix A of the form
A =
1.0000 0.0000 0.0000 0.2345
0.0000 1.0000 0.0000 -0.7653
0.0000 0.0000 1.0000 0.9456
Jacobi Iterative Method
usage: x = jacobi(a, tolerance {optional}, maximum {optional})
jacobi solves the system A*x = b, where A is diagonally dominant.
jacobi will continue until it reaches the maximum number of
iterates or the norm(x(k+1)-x(k)) < tolerance.
INPUT:
a - is an augmented matrix where the coefficients
form a
diagonally dominant matrix,
i.e. |a(i,i)|>sum(a(i,j))
where i <> j.
tolerance - jacobi will continue until the
norm(x(k+1)-x(k)) < tolerance.
Default:
tolerance = max(m,n)*eps*norm(am,"inf").
maximum - the total number of iterations that jacobi will perform
until it quits. Default: maximum
= 100.
OUTPUT:
x - is a matrix whose column vectors are the
solutions to the
system A*x = b
Gauss-Seidel Iterative Method
usage: x = seidel(a, tolerance {optional}, maximum {optional})
seidel solves the system A*x = b, where A is diagonally dominant.
seidel will continue until it reaches the maximum number of
iterates or the norm(x(k+1)-x(k)) < tolerance.
INPUT:
a - is an augmented matrix where the coefficients
form a
diagonally dominant matrix,
i.e. |a(i,i)|>sum(a(i,j))
where i <> j.
tolerance - seidel will continue until the
norm(x(k+1)-x(k)) < tolerance.
Default:
tolerance = max(m,n)*eps*norm(am,"inf").
maximum - the total number of iterations that jacobi will perform
until it quits. Default: maximum
= 100.
OUTPUT:
x - is a matrix whose column vectors are the
solutions to the
system A*x = b.
Forward Gaussian Elimination
forward performs gaussian elimination on a lower triangular matrix.
The matrix should be of size m X m+1.
Usage: x=forward(am)
INPUT:
am - an augmented matrix that should be lower triangular.
EXAMPLE:
a = [2 0 0 -1
1 2 0 5
-2
1 5 6]
OUTPUT:
x - is a column vector of the solutions to the
system.
Backward Gaussian Elimination
back performs gaussian elimination on an upper triangular matrix.
The matrix should be of size m X m+1.
Usage: x=back(am)
INPUT:
am - an augmented matrix that should be upper triangular.
EXAMPLE:
a = [2 1 5 -1
0 2 3 5
0 0 4 6]
OUTPUT:
x - is a column vector of the solutions to the
system.
Tridiagonal matrix reduction.
tridia solves a matrix equation of the form A*x = b where A is a
triadiagonal matrix. A m x 4 matrix is entered where the first
three columns are the diagonals and the last column is b. The
column vector x is returned.
usage: x = tridia(a)
INPUT:
a - is a m x 4 matrix where the first three columns
are the
diagonals of the coefficient
matrix and the last column is
the constant vector, b. Note
that the first element of the
first column is zero and the
last element of the third
column is also zero.
EXAMPLE:
a = [0 1 2 -1
1 2 3 5
-2
1 5 6
3 5 4 -2
6 7 0 5]
OUTPUT:
x - is a column vector of the solutions
to the system.
Octave has a built-in function roots that will find the roots of
a
polynomial equation.
usage: rts = roots(polynomial)
roots will approximate all n roots of a nth degree polynomial equation.
The polynomial must be defined as a row vector of real coefficients.
The
polynomial is of the form a(0)*x^n + a(1)x^(n-1) + a(2)x^(n-2)
+ ... + a(n).
INPUT:
Polynomial - a row vector containing the real coefficients of a
polynomial,
i.e.[a(0), a(1), a(2), ... ,a(n)]. Gaps must be filled with zeros.
OUTPUT:
rts - the approximate roots
for the polynomial equations.
Lagrange Interpolation Method
usage: y = lagrange(x, xdata, ydata)
lagrange will produce interpolated value(s) for an
interpolant or vector of interpolants. The interpolated value(s)
will match p(x) from a polynomial of degree n-1, where n
is the
number of data points.
INPUT:
x - the interpolant value(s) expressed
as a vector.
xdata - a vector of the abscissa of the fixed points.
ydata - a vector of the ordinates of the fixed points.
OUTPUT:
y - a vector of interpolate value(s).
Divided Difference Table
usage: Table=divdiff(xdata,ydata)
divdiff will produce an upper triangular matrix of Newton divided
differences for a set of data points described by two vectors,
xdata and ydata.
INPUT:
xdata - a vector of the abscissa of fixed data points.
ydata - a vector of the ordinates of fixed data points.
OUTPUT:
Table - an upper triangular matrix of Newton divided differences.
Newton Forward Interpolation
usage: [y, c] = newtondivdif(x,xdata,ydata)
newtondivdif will produce interpolated value(s) for an interpolant
vector. The interpolated value(s) will match p(x) from
a polynomial of degree n-1, where n is the number of data points.
The coefficients of the polynomial are derived from Newton divided
differences in a forward direction.
INPUT:
x - the interpolant value(s) expressed
as a vector.
xdata - a vector of the abscissa of fixed data points.
ydata - a vector of the ordinates of fixed data points.
OUTPUT:
y - a vector of interpolate value(s).
c - the coefficients of the interpolating polynomial.
Bezier parametric curve fitting.
usage: [xx yy]=bezier(x,y)
bezier fits a curve to four control points. The curve is contained
in the convex hull if one is defined by the four points. The curve
coincides with the end points.
INPUT:
x - a vector of the abscissa of fixed data points.
y - a vector of the ordinates of fixed data points.
OUTPUT:
xx - the interpolated x values of the curve.
yy - the interpolated y values of the curve.
Natural Spline Interpolation
usage: y = natspline(x, xdata, ydata)
natspline produces interpolated value(s) for an interpolant
vector. The interpolated value(s) will match a cubic polynomial
between the knots, data points. The end knots are left unclamped.
Thus, the slopes at the ends are linear.
INPUT:
x - the interpolant value(s) expressed as a vector.
xdata - a vector of the abscissa of fixed data points.
ydata - a vector of the ordinates of fixed data points.
OUTPUT:
y - a vector of interpolate value(s).
Octave has a built-in function polyfit that finds the best-fit
polynomial for a set of ordered pairs. The best-fit is determined
through regression, the sum of the squared error is minimized.
usage: p=polyfit(x, y, degree)
polyfit will find the polynomial of the degree given, that best
fits the data described by vector x and y. The polynomial
is
described by a vector of its coefficients.
INPUT:
x - a vector of independent data.
y - a vector of dependent data that corresponds
to x.
degree - the degree of the regression polynomial. The degree
should be no more than
n-1, where n is the number of
data points.
OUTPUT:
p - a vector of coefficients that corresponds
to the
regression polynomial.
The coefficients are in
descending order.
EXAMPLE:
p =
1 -2 4 -12
Therefore p(x)= 1x^3 - 2x^2 + 4x - 12.
Octave has a built-in function polyval that evalutes a
polynomial for a value or values.
usage: f=polyval(p, x)
polyval will evaluate a polynomial describe by p at the value(s)
in the vector x.
INPUT:
p - a vector of coefficients for a polynomial.
The
polynomial must be in
descending order.
x - a vector of values at which the polynomial
function
is to be evaluated.
OUTPUT:
f - a vector polynomial function values.
dy
-- = f(t, y).
dt
usage: dfield(diffeq, t0, y0 [optional], t [optional], y [optional])
dfield computes the slope values at integral lattice points and
plots tick
marks corresponding to the slopes of a user supplied function.
The function
contains a single first order ODE.
INPUT:
diffeq - a string name of a m-file containing
a single ODE. The file
must be of the form dydt = fcn (y, t). y is the solution
vector
and t is a scalar (time). dydt is the returned derivative
vector.
EXAMPLE:
function dydt = funn(t, y)
dydt = t.*y.^2;
endfunction
t0 - the initial value
of the independent variable, t.
y0 - the initial condition
of the dependent variable. If omitted,
y0 = t0 and t0 = 0. [optional]
t - [t_min,
t_max]. Range of the independent variable. Default
is [-10, 10]. [optional]
y - [y_min, y_max].
Range of the dependent variable. Default
is [-10, 10]. [optional]
OUTPUT: A plot of a direction field corresponding
to the ODE supplied
by diffeq. An approximate solution is also plotted as long as
more than one input argument is provided.
HINT: Additional solution curves can be plotted by holding
the direction field plot and using diffsolve.
EXAMPLE:
octave: 12> hold on
octave: 13> diffsolve("funn",t0, y0, [-10 10]);
The function euler can be used to solve ODES of the form
dy
-- = f(t, y).
dt
usage: [t_out, y_out] = euler(diffeq, y0, t0, tFinal, step [optional])
euler implements Euler's numerical method of integrating a single
first order
ODE or a system of ODEs.
INPUT:
diffeq - a string name of a m-file containing
a single ODE or system of
ODEs. The file must be of the form dydt = fcn (y, t). y is
the
solution vector and t is a scalar (time). dydt is the returned
derivative vector.
EXAMPLE:
function dydt = funn(t, y)
dydt = t.*y.^2;
endfunction
y0 - the initial condition
of the dependent variable.
t0 - the initial value
of the independent variable, t.
tFinal - the final value of the independent variable,
t.
step - the stepsize of the independent
variable, t. Default value is
(tFinal - t0)/100. [optional]
OUTPUT:
t_out - a column vector containing the
values of the independent variable.
y_out - the numerical approximations for
y(t(i)).
The first row of the
matrix corresponds to y0.
If no output arguments are provided then a graph of the approximate
solution is plotted.
The function heun can be used to solve ODES of the form
dy
-- = f(t, y).
dt
usage: [t_out, y_out] = heun(diffeq, y0, t0, tFinal, step [optional])
heun implements Heun's numerical method of integrating a single
first order
ODE or a system of ODEs.
INPUT:
diffeq - a string name of a m-file containing
a single ODE or system of
ODEs. The file must be of the form dydt = fcn (y, t). y is
the
solution vector and t is a scalar (time). dydt is the returned
derivative vector.
EXAMPLE:
function dydt = funn(t, y)
dydt = t.*y.^2;
endfunction
y0 - the initial condition
of the dependent variable.
t0 - the initial value
of the independent variable, t.
tFinal - the final value of the independent variable,
t.
step - the stepsize of the independent
variable, t. Default value is
(tFinal - t0)/100. [optional]
OUTPUT:
t_out - a column vector containing the
values of the independent variable.
y_out - the numerical approximations for
y(t(i)).
The first row of the
matrix corresponds to y0.
If no output arguments are provided then a graph of the approximate
solution is plotted.
The function midpt can be used to solve ODES of the form
dy
-- = f(t, y).
dt
usage: [t_out, y_out] = midpt(diffeq, y0, t0, tFinal, step [optional])
midpt implements the midpoint numerical method of integrating a
single
first order ODE or a system of ODEs.
INPUT:
diffeq - a string name of a m-file containing
a single ODE or system of
ODEs. The file must be of the form dydt = fcn (y, t). y is
the
solution vector and t is a scalar (time). dydt is the returned
derivative vector.
EXAMPLE:
function dydt = funn(t, y)
dydt = t.*y.^2;
endfunction
y0 - the initial condition
of the dependent variable.
t0 - the initial value
of the independent variable, t.
tFinal - the final value of the independent variable,
t.
step - the stepsize of the independent
variable, t. Default value is
(tFinal - t0)/100. [optional]
OUTPUT:
t_out - a column vector containing the
values of the independent variable.
y_out - the numerical approximations for
y(t(i)).
The first row of the
matrix corresponds to y0.
If no output arguments are provided then a graph of the approximate
solution is plotted.
The function rk4 can be used to solve ODES of the form
dy
-- = f(t, y).
dt
usage: [t_out, y_out] = rk4(diffeq, y0, t0, tFinal, step [optional])
rk4 implements a fourth order Runge-Kutta numerical method of integrating
a single first order ODE or a system of ODEs.
INPUT:
diffeq - a string name of a m-file containing
a single ODE or system of
ODEs. The file must be of the form dydt = fcn (y, t). y is
the
solution vector and t is a scalar (time). dydt is the returned
derivative vector.
EXAMPLE:
function dydt = funn(t, y)
dydt = t.*y.^2;
endfunction
y0 - the initial condition
of the dependent variable.
t0 - the initial value
of the independent variable, t.
tFinal - the final value of the independent variable,
t.
step - the stepsize of the independent
variable, t. Default value is
(tFinal - t0)/100. [optional]
OUTPUT:
t_out - a column vector containing the
values of the independent variable.
y_out - the numerical approximations for
y(t(i)).
The first row of the
matrix corresponds to y0.
If no output arguments are provided then a graph of the approximate
solution is plotted.