interactive numerical computations with octave

octave> version
ans = 4.0.1

shell example

$> octave -q <<NB
> eval ("2 + 4.1")
> NB
ans =  6.1000

в оболочке интерпретатора Octave пользователь может вводить команды языка Octave. если строка заканчивается символом ; то результаты на терминал не выводятся. если в конце строки символ ; отсутствует, то результаты работы выводятся на терминал. текст, который идет после символа % является строкой комментария и интерпретатором не обрабатывается. в системной переменной ans хранится результат последней операции, если команда не содержит знака присваивания. значение переменной ans изменяется после каждого вызова команды без операции присваивания

для обозначения мнимой единицы зарезервировано два имени i, j , поэтому ввод комплексного числа производится в формате: действительная часть + i * мнимая часть или действительная часть + j * мнимая часть

real(Z) выдает действительную часть комплексного аргумента Z
imag(Z) выдает мнимую часть комплексного аргумента Z
angle(Z) вычисляет значение аргумента комплексного числа Z в радианах от –π до π
conj(Z) выдает число комплексно сопряженное Z

< меньше
> больше
== равно
~= не равно
<= меньше или равно
>= больше или равно

octave knows about

concatenation of operations is "."

during the session you can save vars to file (save) and then you can load them from this file (load) S

octave:> load "filenameIn.m" 
octave:> t = 1:0.001:2;
octave:> plot (t, sin(t));
octave:> saveas (1, "filenamePic1", "fmt");
octave:> function f (a)
> 2 * a + 1
> end
octave:> x = [1,2,3,4,5,6];
octave:> y = [0.16, 0.16, 0.16, 0.16,0.16, 0.16];
octave:> bar (x, y);
octave:> saveas (2, "filenamePic2", "fmt");
octave:> save "filenameOut.m"
octave:> quit

the base of the natural logarithm (e) and the ratio of a circle's circumference to its diameter (pi) are pre-defined

to display the value of a variable, simply type the name of the variable at the prompt

you can recall, edit, and reissue previous commands using Emacs/vi style editing commands

once you know the name of the function you wish to use, you can get more help on the function by simply including the name as an argument to help. for example,

    help plot 
will display the help text for the plot function

ending a command with a semicolon tells Octave not to print the result of the command


eps [value]
machine precision

format [ bank | long | short ]
number representaion

clear [variable]

diary on | off

eval (string)
feval ("func_name", args)

terminal input/output

input (prompt)
disp (var)

menu (title, opt1, . . . )
print a title string followed by a series of options
each option will be printed along with a number
the return value is the number of the option selected by the user

a = menu ("example", "say bye", "say good", "quit");

if (a == 3) disp ("quit") 
elseif (a == 1) disp ("bye") 
else disp ("good") 

yes_or_no (prompt)
ask the user a yes-or-no question
return 1 if the answer is 'yes'
takes one argument, which is the string to display
it should end in a space
'yes-or-no-p' adds '(yes or no) ' to it
the user must confirm the answer with RET and
can edit it until it has been confirmed.

just as in C-lang
almost as in C-lang, but instead address length is used as second param

file input/output

save [-ascii | -binary] [variable] filename
save octave data into file
load filename
load octave data from the file

fopen (filename)
returns fid
fscanf (fid, format, count)
read data from file specified by the file descriptor. count may be: inf, NR, [NR,NC], [NR, inf]
fread (fid)
read binary data from the file specified by the file descriptor
fprintf (fid, format, ...)
write data to the file specified by the file descriptor
fclose (fid)
close file with the fid

for example,

myfile = fopen ("splat.dat", "r");

the possible values `mode' may have are
  `r'    for reading
  `w'    for writing. the previous contents are discared
  `a'    open or create a file for writing at the end
  `r+'   open an existing file for reading and writing
  `w+'   open a file for reading or writing. the previous contents are discared
  `a+'   open or create a file for reading or writing at the end

$> cat a.dat 
1 2 3
4 5 6
7 8 9
$> octave -q
octave:1> a = load("a.dat")
a =

   1   2   3
   4   5   6
   7   8   9

octave:2> a
a =

   1   2   3
   4   5   6
   7   8   9

if there's a text header, the file won't load. if you have more than two fields at any point in the file, it won't load. the input file has to have the same number of columns (and each column must have the same number of rows), otherwise octave will complain about the data

> fid = fopen ('data.dat', 'r') ;
> data = fscanf (fid, '%f') ;
> nRows = data (1) ;
> data = reshape (data(2:end), nRows, nRows) ;

you can use the function dlmread():

data = dlmread(file, sep, r0, c0)
which read the matrix data from a text file which uses the delimiter sep between data values
if sep is not defined the separator between fields is determined from the file itself
given two scalar arguments r0 and c0, these define the starting row and column of the data to be read
these values are indexed from zero, such that the first row corresponds to an index of zero

so you can simply use the one-liner:

  data = dlmread ('input_file', ' ', 1, 0);
by calling the function with r0 set to 1, you're effectively skipping the first line, which contains the (now useless) number of rows


to plot line from (a1,b1) to (a2,b2) use

octave:> plot ([3,4],[7,8]) 

octave:> t = linspace(-1,2,100);
octave:> plot (t, arrayfun(@ (x) -.5*(x^2 - x) , t))

to plot a vertical line

octave:> plot ([x; x], [0; 5]) 

to plot a horizontal line

octave:> plot ([0; 40], [y; y]) 

hist (x)
bar (x, y)
stairs (x)
plot (x)
plot (x, y)
plot3 (x, y, z)
color: 'k' blacK 'r' Red 'g' Green 'b' Blue 'm' Magenta 'c' Cyan 'w' White
linestyle: '-' solid lines (default). '--' dased lines. ':' dotted lines. '-.' dash-dotted lines
markerstyle: '+' crosshair 'o' circle '*' star '.' point 'x' cross 's' square 'd' diamond

grid (on/off)
hold (on/off)
subplot(Rows, Cols, Index)

stem (x)
stem (x, y)
polar(θ, ρ)
rose (x)
rose (x,N)
rose (x,y)
contour (Z)
contour (X, Y, Z)
contour3 (Z)
contour3 (X, Y, Z)
mesh (X, Y, Z)
meshc (X, Y, Z)
meshz (X, Y, Z)
surf (X, Y, Z)
surfc (X, Y, Z)
surfc (X, Y, Z, colorbar, shading interp)
waterfall (X, Y, Z)

loglog (x, y)
semilogx (x, y)
semilogy (x, y)

axis ([x_lo, x_hi])
axis ([x_lo, x_hi, y_lo, y_hi])

title (string)
xlabel (string)
ylabel (string)
zlabel (string)
text (X, Y, string)
gtext (string)
legend (string1, string2, ...)

text (10, 80, strcat ("slope is ", num2str(a), "\nintersect is ", num2str(b)));


linspace (begin, end, num)
logspace (begin, end, num)


round (x)
floor (x)
ceil (x)
fix (x)

factorial (n)

gcd (a, b)
lcm (a, b)

exp (x)
log (x)

sqrt (x)
cbrt (x)

complex numbers

abs (z)
arg (z)
angle (z)
real (z)
imag (z)
conj (z)
[theta,r] = cart2pol (x,y)
[x,y] = pol2cart (theta,r)


sin (x)
cos (x)
tan (z)
asin (x)
acos (x)
atan (z)


sinh (x)
cosh (x)
tanh (x)


you can define polynomial as:

> p = [1 2 -1]
which corresponds to x² + 2*x - 1

polyval (p, x)
evaluate the polynomial p at the specified point x

> p = [1, -1, 1];
> polyval (p, 4)
ans =  13

roots (p)
for a vector p with n components, return the roots of the polynomial
p(1) * x^(n-1) + ... + p(n-1) * x + p(n)

> p = [1, -1, 1];
> roots (p)
ans =

   0.50000 + 0.86603i
   0.50000 - 0.86603i

conv (p, q)
the output convolution is a vector with length equal to 'length (p) + length (q) - 1'
when p and q are the coefficient vectors of two polynomials, the convolution represents the coefficient vector of the product polynomial
> p
p =

   2   1

> q
q =

   1   2

> % (2*x + 1) * (x + 2) = 2*x^2 + 5*x + 2
> conv(p, q)
ans =

   2   5   2

[b, r] = deconv (y, a)
solves for b and r such that 'y = conv (a, b) + r

residue (a, b)
[R, P, K, E] = residue (B, A)
[B, A] = residue (R, P, K)
[B, A] = residue (R, P, K, E)

the first calling form computes the partial fraction expansion for the quotient of the polynomials, B and A
K - vector is a polynomial of order N-1 representing the direct contribution,
E -vector specifies the multiplicity of the m-th residue's pole.

> b = [1, 1, 1] ;
> a = [1, -5, 8, -4] ;
> [r, p, k, m] = residue (b, a) ;
> r
r =


> p
p =


> k
k = [](0x0)
> m
m =


which represents the following partial fraction expansion

        s^2 + s + 1       -2        7        3
   ------------------- = ----- + ------- + -----
   s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)

построить интерполяционный полином n-й степени в Octave можно функцией polyfit(x,y,k) (в этом случае k=length(x)-1) для вычисления коэффициентов полинома, после чего с помощью функции polyval(A,t) вычислить значение полинома в необходимых точках


в Octave есть встроенные функции вычисления суммы (sum), произведения (prod) элементов массива, поиска максимума (max) и минимума (min), сортировки (sort)

length (v)

dot (v, u)
cross (v, u)

norm (v)
norm (v, m)
m = 1 | 2 (default) | inf

diff (v)
grad (v)


.+   .-   .*   ./   .^
element-wise summation, subtruction, multiplication, division, power

sum (v)
prod (v)

> f = [3:7]'
f =


vector fields

gradient (FX,FY,FZ)
divergence (FX,FY,FZ)
curl (FX,FY,FZ)


значение нижней границы индексации в Octave равно единице

ввод элементов матрицы осуществляется в квадратных скобках, при этом элементы строки отделяются друг от друга пробелом или запятой, а строки разделяются между собой точкой с запятой. обратиться к элементу матрицы можно, указав после имени матрицы, в круглых скобках, через запятую, номер строки и номер столбца на пересечении которых элемент расположен

to create a new matrix and store it in a variable:

octave:> A = [ 1, 1, 2; 3, 5, 8; 13, 21, 34 ] ;
Octave uses a comma or space to separate entries in a row, and a semicolon or carriage return to separate one row from the next. for example,
      B = rand (3, 2); 
will create a 3 row, 2 column matrix with each element set to a uniformly distributed random value between zero and one
length (A)
returns maximal dimension
size (A)
returns all dimensions
columns (A)
rows (A)

det (A)
trace (A)
rank (A)


zeros (n)
returns square matrix with all els = 0
zeros (i, j)
returns matrix i x j with all els = 0
ones (n)
returns square matrix with all els = 1
ones (i, j)
returns matrix i x j with all els = 1
eye (n)
returns nxn identical matrix
eye (i, j)
returns matrix i x j with all main diag elems = 1

rand (n, m)
generate matrix nxm with each elemet ~ U[0,1]
randn (n, m)
generate matrix nxm with each element ~ N(0,1)
rand (k)
generate kxk matrix with each element ~ U[0,1]
randn (k)
generate kxk matrix with each element ~ N(0,1)

blkdiag (A, B, C, ...)
build a block diagonal matrix from A, B, C, ...
all arguments must be numeric and either two-dimensional matrices or scalars


vertcat (A,B)
return the vertical concatenation of array objects
horzcat (A,B)
return the horizontal concatenation of array objects


fliplr (x)
return a copy of x with the order of the columns reversed
flipud (x)
return a copy of x with the order of the rows reversed


A (i,:)
take i row of A

A (:,j)
take j column of A

A ([2 3],:) = A([3 2], :)
exchanges columns 2 and 3 of A

A (:, [2 3]) = A(:, [3 2])
exchanges rows 2 and 3 of A

A (2:4, 3:7)
returns rows 2 to 4 and cols 3 to 7 (3x5 matrix)


.+   .-   .*   ./   .^
element-wise summation, subtruction, multiplication, division, power

A > N
returns matrix of booleans
A < N
returns matrix of booleans

functions like exp, log, sin, cos, sqrt, etcetera, can also be applied to matrices. in this case, they are simply applied to every element of the matrix

the first and the second also have a linear algebraic counterpart: the matrix exponential and matrix logarithm. the names of the corresponding functions is formed with the suffix m, so the matrix exponential is denoted by expm, matrix logarithm by logm

exp (A)
result is the matrix B with elements: bi,j = e^ai,j
expm (A)
result is the matrix B: B = I + A + A*A ./ 2! + A*A*A ./ 3! + ...
A \ b
if the system A * x = b has a unique solution, then A\b returns this solution
if there is no solution, then A\b computes the least-square solution
if there is more than one solution, then A\b just returns an arbitrary solution
b \ A
returns x: x * A = b
kron (A, B)
tensor product of two matrices

transpose the matrix
inv (A)
returns inverse
pinv (A)
returns pseudo-inverse

null (A)
returns a basis for the nullspace
orth (A)
computes an orthonormal basis for the range

norm (A)
computes the L₂ norm
norm (A, m)
m = 1 | 2 (default) | inf

poly (A)
computes the characteristic polynomial det (λ * I − A)
eigenvalues of A correspond to zeroes of this polynomial
cond (A)
the condition number is the factor, with which perturbation in the input are magnified:
condition number = error in output/error in input
if the matrix A is orthogonal, then the condition number of A is 1


[S, E] = eig (A)
gives two matrices: S - eig-vects for matrix A, E - diag - eig-vals for matrix A

[L, U, P] = lu (A)
P * A = L * U, L - low-triangular matrix, U - upper triangular matrix, P - permutation matrix

[U, S, V] = svd (A)
[U, S, V] = svd (A, ECON)
it factorizes a possibly nonsquare matrix A in three matrices, such that U * S * V’ = A
one of the reasons why the singular value decomposition is in wide use: it decomposes a matrix A in the product A = U * Σ * V⁻ , where the matrices U and V are orthogonal. this is leads often to better algorithms than the eigendecomposition A = V * Λ * V⁻, where the matrix V is not orthogonal

[Q, R] = qr (A)
gives orthogonal Q and triangular R, Anm = Qnn * Rnm

C = chol (A)
Cholesky decomposition symmetric positive definite matrix A: A = C⁺ * C
btw - it is the best way to know is A is positive definite or not

[P, H] = hess (A)
compute the Hessenberg decomposition of the matrix A
P * H * P = A where
P is a square unitary matrix: P' * P = I
H(i,j) = 0 ∀ i > = j+1

something else

tril (A)
return a new matrix formed by extracting the lower triangular part of the matrix A, and setting all other elements to zero
triu (A)
return a new matrix formed by extracting the upper triangular part of the matrix A, and setting all other elements to zero

diag (v)
returns matrix with diagonal as vector v
diag (A)
retruns vector consists from diag elems of matrix A

toeplitz (c)
toeplitz (c, r)
return the Toeplitz matrix constructed from the first column C, and optionally the first row R
vander (c, n)
returns the Vandermonde matrix whose next to last column is c if n is specified, it determines the number of columns; otherwise, n is taken to be equal to the length of c


union (a, b)
intersect (a, b)
ismember (x, a)
returns 0 | 1


min (x)
max (x)

mean (x)
if x is a vector, compute the mean of the elements of x
if x is a matrix, compute the mean for each column and return them in a row vector
median (x)
if x is a vector, compute the median value of the elements of x
if x is a matrix, compute the median value for each column and return them in a rowvector
var (x)
for vector arguments, return the (real) variance of the values
for matrix return variance for each column
std (x)
if x is a vector, compute the standard deviation of the elements of x
if x is a matrix, compute the standard deviation for each column

cov (x, y)
if each row of x and y is an observation and each column is a variable,
the (i, j)-th entry of cov (x,y) is the covariance between the i-th variable in x and the j-th variable in y
corr (x, y)
the (i,j)-th entry of cor(x,y) is the correlation between the i-th variable in x and the j-th variable in y for matrices


[v, ier, nfun, err] = quad (f, a, b, tol, sing)
numerical integration based on Gaussian quadrature

f - the name of the function, the function handle or the inline function.
a and b - limits of integration. either or both may be infinite
the optional argument tol is a vector that specifies the desired accuracy of the result. the first element of the vector is the desired absolute tolerance, and the second element is the desired relative tolerance
the optional argument sing is a vector of values at which the integranl is known to be singular

the result of the integration is returned in v
ier contains an integer error code (0 indicates a successful integration)
nfun indicates how many function evaluations were required
err contains an estimate of the error in the solution

имя функции должно быть в кавычках, потому что этот параметр - строка


вычислим длину дуги синусоиды на интервале полного периода и сравним с отрезком прямой [0,2p]

в файле arc.m

function y = arc (x)
  y = sqrt (1 + sin (x) * sin (x)) ;

теперь в репле октава :

octave> z = quad ("arc", 0, 6.2832) ;
octave> z / 6.2832 
ans =  1.2160  

double integrals

dblquad (f, xa, xb, ya, yb, tol, quadf, . . . )
numerically evaluate a double integral

the function over with to integrate is defined by f, and the interval for the integration is defined by [xa, xb, ya, yb]. the function f must accept a vector x and a scalar y, and return a vector of the same length as x. If defined, tol defines the absolute tolerance to which to which to integrate each sub-integral. additional arguments, are passed directly to f

to use the default value for tol one may pass an empty matrix

ordinary differential equations

              dx/dt = f(x, t)
              x(0) = x₀

octave has built-in functions for solving nonlinear differential equations of the form dx/dt = f (x, t) with the initial condition x(t0) = x0 . you must provide a definition of the function f (x, t). this may be accomplished by entering the function body directly on the command line:

octave:> function [y] = f (x, t)
>  r = 0.25;   k = 1.4;
>  a = 1.5;    b = 0.16;
>  c = 0.9;    d = 0.8;
>  y(1) = r * x(1) * (1 - x(1) / k) - a * x(1) * x(2) / (1 + b * x(1));
>  y(2) = c * a * x(1) * x(2) / (1 + b * x(1)) - d * x(2);
> endfunction

> x0 = [1; 2];
> t = linspace (0, 50, 200);
> x = lsode ("f", x0, t);
> plot (t, x);
> print -deps foo.eps
it will create a file called foo.eps that contains a rendering of the current plot in Encapsulated PostScript format

octave can handle both 'stiff' and 'nonstiff' equations

lsode expects three arguments:

the function returns a matrix of x in times, specified by user
options can be specified by lsode_options function

[A, istate, msg] = lsode (fcn, x0, t0, t_crit)
the first argument, fcn, is a string that names the function f to call to compute the vector of rhs for the set of equations
the function f must have the form ydot = f (y, t) in which ydot and y are vectors and t is a scalar
if fcn is a two-element string array or a two-element cell array of strings, inline functions, or function handles, then the first element names the function f described above, and the second element names a function to compute the Jacobian of f. the Jacobian function must have the form jac = j (x, t) in which jac is the matrix of partial derivatives
the second and third arguments specify the initial state of the system, x0 , and the initial value of the independent variable t0. the first element of t should be t0 and should correspond to the initial state of the system x0, so that the first row of the output is x0
the fourth argument is optional, and may be used to specify a set of times that the ODE solver should not integrate past. it is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative

after a computation, each row of the matrix A corresponds to an element of the vector t
after a successful computation, the value of istate will be 2 (consistent with the Fortran version of Lsode). if the computation is not successful, istate will be something other than 2 and msg will contain additional information



% equation of moving in 3Dspace
% dr/dt = v(t) = [x/t, y/t, z/t]


function v = some_func (x, t) 
  d = exp (-t) ;
  v = (eye (3) * d) * x ;

function diff3 (fdot, x0, t)
  f = lsode (fdot, x0, t) ;
  plot (t, f, 'linewidth', 2) ; grid ;
  legend (" x " , " y ", " z ") ;
  xlabel ("time") ; ylabel ("position") ;

x0 = [1; 2; 3] ;
t = linspace (1, 10, 200) ;
diff3 (@some_func, x0, t) ; pause ; 

solution of high-order ODEs

if your problem is of order 2 or higher: rewrite your problem as a first order system. for example : the angle y of an undamped pendulum with a driving force sin(5*t) satisfies the differential equation

    y" = -sin(y) + sin(5*t)   

and the initial conditions:

    y (0) = 1
    y'(0) = 0

now, let y₁=y and y₂=y'

this gives you the first order system:

    y₁' = y₂
    y₂' = -sin(y₁) + sin(5*t)

with the initial conditions

    y₁(0) = 1
    y₂(0) = 0

algebraic differential equations

    f (x, x', t) = 0
    x (0)        = x₀
    x'(0)        = x'₀

the function is called dassl, and it expects four arguments:

the function returns two matrices containing the x and x' at the times specified by the user

[A, B, istate, msg] = dassl (fcn, x0, xdot0, t, t_crit)
the solution is returned in the matrices A and B, with each row in the result matrices corresponding to one of the elements in the vector t

the first element of t should be t0 and correspond to the initial state of the system x0 and its derivative xdot0, so that the first row of the output A is x0 and the first row of the output B is xdot0

the first argument, fcn, is a string, that names the function f to call to compute the vector of residuals for the set of equations. it must have the form res = f (x, xdot, t) in which x, xdot, and res are vectors, and t is a scalar



function yourdif (fdot, x0, d0, t)
  [A, B, c, d] = dassl (fdot, x0, d0, t) ;
  plot ( t , A , "r" , 'linewidth' , 2 ,    t , B , "b" , 'linewidth' , 2 ) ;
  legend ("func" , "deriv") ; pause ;

function r = some_func_here (x, d, t) r(1) = 3 - d ; end

x0 = 5 ; d0 = 3 ; t = linspace (1, 10, 200) ;
yourdif ("some_func_here", x0, d0, t) ;

fast Fourier transforms

classic example: amplitude modulation (AM)

defining a sinusoidal signal:

freq1 = 100;
period1 = 1 / freq1;
w1 = 2 * pi * freq1;
num_tsteps = 1000;
num_periods = 2;
tstep = num_periods * period1 / num_tsteps;
t = 0:tstep:(num_periods * period1);
x1 = sin(w1 * t);
plot(t, x1);

modulate it with a 2.5kHz signal:

freq2 = 2500;
period2 = 1 / freq2;
w2 = 2 * pi * freq2;
x2 = sin(w2 * t);

x3 = x1 .* x2;
plot(t, x3);

the function

fft (a, n)
computes a Fast Fourier Transform (FFT) with n points based on signal samples contained in the list a

the list returned from the function fft isn't dependent on any particular timescale (provided the signal is sampled at a constant rate, called the sampling frequency). it's up to you to set a timescale that makes sense

you know the list from fft has n points - the list of associated frequencies will therefore contain n points. since the signal was sampled with a constant f[s], the spacing between the frequencies will also be equal

now I declare the frequency list as a range of n numbers, which will then be scaled by an appropriate amount:

f = (0:n-1) * scale;

what's the scaling factor? it's easiest to first normalize the frequencies so they run from 0 to 1 (non-inclusive) by dividing by n

the sampling frequency is the maximum frequency because one period of the Fourier transform represents an amount of time equal to the sampling period (1/f_s). I can therefore rewrite the last statement as:

f = f_s * (0 : n-1)/n;

in our case, f_s = num_tsteps / (period1 * num_periods). pulling all of this together, let's define list y as the Fourier transform of the AM signal, and f as its associated frequencies:

n = 2**16;
y = abs (fft (x3, n));

f_s = num_tsteps / (period1 * num_periods);
f = f_s * (0 : n-1) / n;

the list returned from fft contains complex numbers, but we are only interested in the magnitude of each frequency component of the signal

let's plot the FFT:

axis ([0, 5000]);
plot (f, y);

this symmetry, however, is not a reflection across f = 0 Hz, as it should be. to view the true spectrum it's necessary to use the function fftshift (v), which shifts the values in v, a list containing the results of a Fourier transform, such that the frequency 0 Hz is at the center of the list. you'll also have to shift the corresponding frequencies:

y = abs (fft (x3, n));
y = fftshift (y);
f = f_s * (-n/2 : n/2-1) / n;
axis ([-5000, 5000]);
plot (f, y);

there are more fft related functions

ifft (a, n)
converts the inverse Fourier transform (frequency domain -< time domain).

fftconv (a, b, n)
computes the convolution of signals a and b using an n-point FFT (which is significantly faster than conv((a, b)).

fftfilt (b, x, n)
filters signal x using a FIR filter

injecting some random noise by using the randn (rows, columns), which returns a matrix of random numbers with normal distribution

alternatively, we could use rand (rows, columns), but for this purpose we want values centered about zero

let's plot the new modulated signal:

x3 = x3 + randn (1, length(x3));
plot(t, x3);

another examples

%% defining a sinusoidal signal ... 

freq1 = 100 ; period1 = 1 / freq1 ; w1 = 2 * pi * freq1 ;
num_steps = 1000 ; num_periods = 4 ;
dt = num_periods * period1 / num_steps ;
t = 0:dt:dt * num_steps ;
x1 = sin (w1 * t) ; plot (t, x1, 'linewidth', 2) ; pause

%% ... and modulate it with a 250 Hz signal
freq2 = 250 ; period2 = 1 / freq2 ; w2 = 2 * pi * freq2 ;
x2 = sin (w2 * t) ;
x3 = x1 .* x2 ;
plot (t, x3, 'linewidth', 2) ; pause

%% ... and its spectrum now will be 
n = 2^8 ;

y = fft (x3, n) ; df = 1/dt ; f = df *  (0:n-1) / n ;
plot (f, y, 'linewidth', 2) ; pause

LP with packet 'glpk'

octave can be used to solve a linear program (LP)

LP is finding minimax for functions like that:

      f(x₁,x₂,...) = c₁*x₁ + c₂*x₂ + ..

      a₁₁*x₁ + a₁₂*x₂ + ... ≥|≤ b₁
      a₂₁*x₁ + a₂₂*x₂ + ... ≥|≤ b₂

            0  <   xₙ   <  ∞
or      lower  <   xₙ   <  ∞
or      lower  <   xₙ   <  upper
or          0  <   xₙ   <  upper

that's all about it. example:

> A = [2,1;1,1;1,0]
   A =

   2   1
   1   1
   1   0

> c = [3,2]
   c =

   3   2

> b = [100,80,40]
   b =

   100    80    40

> [x_max, z_max] = glpk (c, A, b, [0,0], [], "UUU", "CC", -1)
   x_max =


   z_max =  180

there is glpk function to do the job:

[x_max, z_max] = glpk (c, A, b, l, u, "UUL", "II", -1)
[x_min, z_min] = glpk (c, A, b, l, u, "LUL", "CC", 1)
A - matrix of constraints coeffs
b - vector of constraints rhs for each var
c - vector of objective coeffs for each var
l - vector of low bounds for each var
u - vector of upper bounds for each var
F U S L D - Free(the constraint is ignored), Upper(A*x<b), Eq(A*x=b), Low (A*x>b), Both (bl<A*x<bu)
C I - continious/integer
1 | -1 - minimization|maximization

LP : minimax   c⁺ * x   subject to      a1i * x    cond    b1,
                                        ami * x    cond    bm,

                 x ∈ ℝⁿ
                 c ∈ ℝⁿ
                 b ∈ ℝⁿ

QP with function 'qp'

octave can be used to solve a quadratic program (QP)

function qp solves the quadratic program

[x, OBJ, INFO, Lambda] = qp (x0, H)
[x, OBJ, INFO, Lambda] = qp (x0, H, q)
[x, OBJ, INFO, Lambda] = qp (x0, H, q, A, b)
[x, OBJ, INFO, Lambda] = qp (x0, H, q, A, b, lb, ub)
[x, OBJ, INFO, Lambda] = qp (x0, H, q, A, b, lb, ub, A_LB, A_IN, A_UB)

QP:    minimize       .5  *  x' * H * x  +  x' * q

       subject to      A * x = b
                       lb ≤ x ≤ ub
                       A_lb ≤ A_in * x ≤ A_ub

any bound (A, B, lb, ub, A_LB, A_UB) may be set to the empty matrix [] if not present. if the initial guess is feasible the algorithm is faster

 #!/usr/bin/octave -q
% real integrator with damping

x = zeros (64);
y = zeros (64);

x (2) = 1;
for t = 2:64
  y (t)  =  x (t) + .9 * y (t - 1);

plot (1:64, y, 'linewidth', 2)


функция имеет заголовок вида

    function name1 [, name2, …] = fun(var1 [, var2, …])
где name1 [, name2, …] – список выходных параметров, то есть переменных, которым будет присвоен конечный результат вычислений, fun – имя функции, var1 [, var2, …] – входные параметры. таким образом, простейший заголовок функции выглядит так:
    function name = fun(var)

все имена переменных внутри функции, а так же имена из списка входных и выходных параметров воспринимаются системой как локальные, то есть эти переменные считаются определенными только внутри функции

условный оператор имеет вид:

    if условие

сокращенный условный оператор записывают так:

    if условие

расширенный условный оператор:

    if условие1
    elseif условие2
    elseif условие 3
    elseif условие n

еще одним способом организации разветвлений является оператор альтернативного выбора:

    swith параметр
    case значение1 операторы1
    case значение2 операторы2
    case значение3 операторы3
    otherwise операторы

оператор цикла с предусловием:

    while выражение

для цикла с известным числом повторений:

    for параметр = начальное значение:шаг:конечное значение

операторы передачи управления принудительно изменяют порядок выполнения команд. в языке программирования Octave таких операторов два. операторы break и continue используют только внутри циклов. оператор break осуществляет немедленный выход из циклов while, for и управление передается оператору, находящемуся непосредственно за циклом. оператор continue начинает новую итерацию цикла, даже если предыдущая не была завершена

function Outputs = function_name (Par1, Par2, ...)




function [z] = vect_mul (x , y) 
  d = length (x) ; 
  if (d != length (y))  
    printf ("vectors are not equal\n") ; 
    z = [0] ; 
    return ; 
    for i = 1 : d 
      z(i) = x(i) * y(i) ;

octave> function last_words ()
>         disp ("bye bye kanzas");
>       endfunction
octave> atexit ("last_words");


octave:> f = @(x) x * x + 2 ;
octave:> f(5)
ans = 27

 arrayfun (@ (x) 1/x, Param) 


you can pass a function handle (@function_name) to another function:
function r = computeit (f, g, v)
  r = f (g (v));
computeit (@exp, @sin, pi/3)
computeit (@log, @cos, pi/6)

or pass the string name of the function and use the feval primitive.

function r = computeit2 (f, g, v)
  r = f (feval (g, v));
computeit2(@exp, "sin", pi/3)


arrayfun (FUNC, A)
execute a function on each element of an array
this is useful for functions that do not accept array arguments
the first input argument FUNC can be a string, a function handle, an inline function, or an anonymous function
the second argument A can be a logic array, a numeric array, a string array, a structure array, or a cell array
by a call of the function 'arrayfun' all elements of A are passed on to the named function FUNC individually


a script file is a file containing (almost) any sequence of Octave commands. it is read and evaluated just as if you had typed each command at the Octave prompt, and provides a convenient way to perform a sequence of commands that do not logically belong inside a function it is possible to define more than one function in a single script file and load (but not execute) all of them at once

to do this, the first token in the file (ignoring comments and other white space) must be something other than function. if you have no other statements to evaluate, you can use a statement that has no effect, like this:

# prevent Octave from thinking that this is a function file:


# define function one:

function one ()

to run this script just input filename in octave prompt

if the first token in a file (ignoring comments) is function, Octave will compile the function and try to execute it

unlike a function file, a script file must not begin with the keyword function. if it does, Octave will assume that it is a function file, and that it defines a single function that should be evaluated as soon as it is defined

a script file also differs from a function file in that the variables named in a script file are not local variables, but are in the same scope as the other variables that are visible on the command line

%%   octave --persist lprg_1.m 
%% or
%%   octave:1> eval lprg_1 

t = 0:1:30;

function y = a (t)
  y = 50 - 2 * t;

function y = b (t)
  y = 40 - 4 * t / 3;

function y = c (t)
  y = 50 - 6 * t / 5;

hold; grid; xlabel ("factory A, days"); ylabel ("factory B, days");
plot (t, a(t), 'r', 'linewidth', 2)
plot (t, b(t), 'g', 'linewidth', 2)
plot (t, c(t), 'm', 'linewidth', 2)
plot([24,24],[0,24], 'linewidth', 2) 
plot([0,24],[24,24], 'linewidth', 2) 


the easiest/most_appropriate way to load data from a .mat file onto the workspace is via the load command. it allows you to import only a single variable (whose name you know) into the workspace

you can do this by simply running the load command, without assigning to a variable:

 load ('file1.mat', 'M11');
 load ('file2.mat', 'M21');
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        M11         1x3                         24  double
        M21         1x3                         24  double

Total is 6 elements using 48 bytes
 isequal (M11, M21)
ans =  1

however, if you do collect into a variable, this variable becomes a struct, whose fieldnames correspond to the names of the variables you imported, e.g.

 S1 = load ('file1.mat', 'M11');
 S2 = load ('file2.mat', 'M21');
 isequal (S1.M11, S2.M21)
ans =  1

file octaverc:

   set (0, "defaultlinelinewidth", 2)
save ('my_environment')
load ('my_environment')
saves/restore all the variables in the environment into/from the file
history -w 'my_history'
history -r 'my_history'
command with the -w tag saves the command history to the file


octave> pkg prefix /usr/local/lib/octave/packages/4.0.0 
octave> pkg install some-pkg-name.tar.gz
octave> pkg install --forge some-pkg-name.tar.gz
octave> pkg list
octave> pkg uninstall some-pkg-name