interactive numerical computations with plotting capabilities

``````octave> version
ans = 4.0.1
``````

shell example

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

octave knows about arithmetic operations (+ ,- , *, /), exponentiation (^), natural logarithms/exponents (log, exp), and the trigonometric functions (sin, cos, tan)

octave calculations work on real or imaginary numbers (i,j)

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.

general

eps [value]
machine precision

format [ bank | long | short ]
number representaion

clear [variable]

clc
date
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")
endif;

quit
```

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.

printf
just as in C-lang
scanf
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 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]
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
`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
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

plotting

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)
clf
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)));

intervals

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)

polynomes

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 =

-2.0000
7.0000
3.0000

> p
p =

2.00000
2.00000
1.00000

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

1
2
1
```

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

vectors

length (v)

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

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

diff (v)

map

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

sum (v)
prod (v)

```> f = [3:7]'
f =

3
4
5
6
7
```

vector fields

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

matrices

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)

construction

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

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)
uniform distribution from 0 to 1
randn (n, m)
with mean 0 and sigma 1

A (i,:)
take i row of A

A (:,j)
take j column of A

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 ([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)

by-elem-ops

.+   .-   .*   ./   .^
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'
transpose the matrix
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

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

decompositions

[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)
triu (A)
diag (v)
returns matrix with diagonal as vector v
diag (A)
retruns vector consists from diag elems of matrix A

toeplitz (v)
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

statistics

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

rand (k)
uniform distribution from 0 to 1
randn (k)
normal distribution mean = 0, sigma = 1

functions

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

return

exit

exit(Status)

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

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

lambda

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

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

HOF

you can pass a function handle (@function_name) to another function:
``````function r = computeit (f, g, v)
r = f (g (v));
endfunction

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));
endfunction

computeit2(@exp, "sin", pi/3)
``````

map

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

integrals

[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

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 name of the solver
• the initial conditions
• a vector of timestamps

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

example

``````#!/home/guest/.nix-profile/bin/octave

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

1;

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

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

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 name of the solver
• the initial position
• the initial speed
• a vector of times
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

example

``````#!/home/guest/.nix-profile/bin/octave

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 ;
end

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

```> 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 =

20
60

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);
end;

plot (1:64, y, 'linewidth', 2)
pause;
``````

scripts

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:

1;

# 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;
end

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

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

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

environment

file octaverc:
```   set (0, "defaultlinelinewidth", 2)
```
save ('my_environment')
```octave> pkg prefix /usr/local/lib/octave/packages/4.0.0