\documentclass[a4paper]{article}
\usepackage[latin1]{inputenc}
\usepackage{graphicx,amsmath}

\newcommand{\Octave}{{\textsf{Octave}}}
\newcommand{\Matlab}{{\textsf{Matlab}}}

\begin{document}

\title{Introduction to GNU Octave}
\date{2008/08/16}
\author{Hubert Selhofer, revised by Marcel Oliver \\ updated to
current \Octave\ version by Thomas L. Scofield}
\maketitle

\begin{center}
  \includegraphics[height=10cm]{sombr} %.pdf
\end{center}

\tableofcontents


\section{Basics}

\subsection{What is \Octave?}

\Octave{} is an interactive programming language specifically suited
for vectorizable numerical calculations.  It provides a high level
interface to many standard libraries of numerical mathematics, e.g.\ 
LAPACK or BLAS.

The syntax of \Octave{} resembles that of \Matlab{}.  An \Octave{}
program usually runs unmodified on \Matlab{}.  \Matlab, being
commercial software, has a larger function set, and so the reverse
does not always work, especially when the program makes use of
specialized add-on toolboxes for \Matlab.

\subsection{Help!}

\begin{itemize}
\item \texttt{help} lists the names of all built-in functions and
internal variables.
\item \texttt{help name} further explains the variable or function
``name''.
\end{itemize}

\paragraph*{Example:}
\begin{verbatim}
  octave:1> help eig
\end{verbatim}


\subsection{Input conventions}

\begin{itemize}
\item All commands can be typed in at the prompt or read from a
script.
\item Scripts are plain text files with file suffix \texttt{.m}.  They
are imported by calling the file name without the suffix and behave as
if their content was typed in line by line.
\item \texttt{;} separates several commands within a line.  A command
terminated by \texttt{;} does not display its result on-screen.
\item \texttt{,} separates two commands without suppressing on-screen
output.
\item \texttt{...} at the end of the line denotes that an expression
continues into the next line.
\item Comments are preceded by \texttt{\%}.
\item \Octave{} is case sensitive.
\end{itemize}


\subsection{Variables and standard operations}

\begin{itemize}
\item \texttt{varname = expression} assigns the result of
\texttt{expression} to \texttt{varname}.
\item \Octave\ has all the usual mathematical functions \texttt{+},
\texttt{-}, \texttt{*}, \texttt{/}, \texttt{\^{}}, \texttt{sin},
\texttt{cos}, \texttt{exp}, \texttt{acos}, \texttt{abs}, etc.
\item The operators of elementary logic are:

\begin{tabular}{clclcl}
  \texttt{<} & smaller \qquad &
  \texttt{<=} & smaller or equal \qquad &
  \texttt{\&} & and \\
  \texttt{>} & greater &
  \texttt{>=} & greater or equal &
  \texttt{|} & or \\
  \texttt{==} & equal &
  \verb+~=+ & not equal &
  \verb+~+ & not
\end{tabular}
\end{itemize}
When debugging user-defined objects, the following commands are
useful:
\begin{itemize}
\item \texttt{whos} shows all user-defined variables and functions
(see Section~\ref{functions}).
\item \texttt{clear name} clears \texttt{name} from memory; if no name
is given, \emph{all} variables and functions will be cleared.
\item \texttt{type name} displays information about the variable or
function \texttt{name} on the screen.
\end{itemize}

\paragraph*{Examples:}
\begin{verbatim}
  octave:1> x12 = 1/8, long_name = 'A String'
  x12 = 0.12500
  long_name = A String
  octave:2> sqrt(-1)-i
  ans = 0
  octave:3> x = sqrt(2); sin(x)/x
  ans = 0.69846
\end{verbatim}
And here is a script \texttt{doless}, saved in a file named
\texttt{doless.m}:
\begin{verbatim}
  one = 1;
  two = 2;
  three = one + two;
\end{verbatim}
Calling the script:
\begin{verbatim}
  octave:1> doless
  octave:2> whos
  *** local user variables:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   real scalar           1      1  three
   wd   real scalar           1      1  one
   wd   real scalar           1      1  two
\end{verbatim}


\section{Vector and matrix operations}

Matrices and vectors are the most important building blocks for
programming in \Octave.

\subsection{Vectors}
\label{vector}

\begin{itemize}

\item Row vector: $v = (1,2,3)$:

  \texttt{v = [ 1 2 3 ]}

\item Column vector: $v = (1,2,3)^T$:

  \texttt{v = [ 1; 2; 3 ]}

\item Automatic generation of vectors with constant increment:

  \texttt{Start[:Increment]:End}
\end{itemize}

\paragraph*{Examples:}
\begin{verbatim}
  octave:1> x = 3:6
  x =
    3  4  5  6
  octave:2> y = 0:.15:.7
  y =
  0.00000  0.15000  0.30000  0.45000  0.60000
  octave:3> z = pi:-pi/4:0
  z =
  3.14159  2.35619  1.57080  0.78540  0.00000
\end{verbatim}


\subsection{Matrices}

A matrix $\displaystyle A = \begin{pmatrix} 1 & 2 \\ 3 & 4
\end{pmatrix}$ is generated as follows.
\begin{verbatim}
  octave:1> A = [ 1 2; 3 4]
  A =
    1  2
    3  4
\end{verbatim}
Matrices can assembled from submatrices:
\begin{verbatim}
  octave:2> b = [5; 6];
  octave:3> M = [A b]
  M =

    1  2  5
    3  4  6
\end{verbatim}

There are functions to create frequently used $m \times n$ matrices.
If $m=n$, only one argument is necessary.
\begin{itemize}
\item \texttt{eye(m,n)} produces a matrix with ones on the main
diagonal and zeros elsewhere.  When $m=n$, the identity matrix is
generated.
\item \texttt{zeros(m,n)} generates the zero matrix of dimension $m
\times n$.
\item \texttt{ones(m,n)} generates an $m\times n$ matrix where all
entries are $1$.
\item \texttt{rand(m,n)} generates a random matrix whose entries are
uniformly distributed in the interval $(0,1)$.  
\end{itemize}


\subsection{Basic matrix arithmetic}

\begin{itemize}
\item \texttt{+}, \texttt{-}, and \texttt{*} denote matrix addition,
subtraction, and multiplication.
\item \texttt{A'} transposes and
conjugates \texttt{A}.
\item \texttt{A.'} transposes \texttt{A}.
\end{itemize}

\paragraph*{Examples:}
\begin{verbatim}
  octave:1> A = [1 2; 3 4]; B = 2*ones(2,2);
  octave:2> A+B, A-B, A*B
  ans =
    3  4
    5  6
  ans =
    -1   0
     1   2
  ans =
     6   6
    14  14
\end{verbatim}


\subsection{Element-wise operations}
\label{element}

While \texttt{*} refers to the usual matrix multiplication,
\texttt{.*} denotes element-wise multiplication.  Similarly,
\texttt{./} and \texttt{.\^{}} denote the element-wise division and
power operators.

\paragraph*{Examples:}
\begin{verbatim}
  octave:1> A = [1 2; 3 4]; A.^2   % Element-wise power
  ans =

     1   4
     9  16

  octave:2> A^2   % Proper matrix power: A^2 = A*A
  ans =

     7  10
    15  22
\end{verbatim}


\subsection{Indexing and slicing}

\begin{itemize}
\item \texttt{v(k)} is the $k$-th element of the row or column vector
$v$.
\item \texttt{A(k,l)} is the matrix element $a_{kl}$.
\item \texttt{v(m:n)} is the ``slice'' of the vector $v$ from its
$m$th to its $n$-th entry.
\item \texttt{A(k,:)} is the $k$th row of matrix $A$. \item
\texttt{A(:,l)} is the $l$th column of matrix $A$.
\end{itemize}
Querying dimensions:
\begin{itemize}
\item \texttt{length(v)} returns the number of elements of the vector
$v$.
\item \texttt{[Rows,Columns] = size(A)} returns the number of rows and
columns of the matrix $A$.
\end{itemize}
Reshaping:
\begin{itemize}
\item \texttt{reshape(A,m,n)} transforms \texttt{A} into an $m \times
n$-matrix.
\item \texttt{diag(A)} creates a column vector containing the diagonal
elements $a_{jj}$ of the matrix \texttt{A}.
\item \texttt{diag(v)} generates a matrix with the elements from the
vector \texttt{v} on the diagonal. 
\item \texttt{diag(v,k)} generates a matrix with the elements from the
vector \texttt{v} on the \texttt{k}-th diagonal.
\item \texttt{A(k,:)\ = rv} assigns to the \texttt{k}-th row of
\texttt{A} the row vector \texttt{rv}.
\item \texttt{A(k,:)\ = []} deletes the \texttt{k}-th row of
\texttt{A}.
\item \texttt{A(:,j) = cv} assigns to the \texttt{j}-th column of
\texttt{A} the column vector \texttt{cv}.
\item \texttt{A(:,j) = []} deletes the \texttt{j}-th column of
\texttt{A}.
\end{itemize}

\paragraph*{Examples:}
\begin{verbatim}
  octave:1> A = [1 2 3; 4 5 6]; v = [7; 8];
  octave:2> A(2,3) = v(2)
  A =

    1  2  3
    4  5  8

  octave:3> A(:,2) = v
  A =

    1  7  3
    4  8  8

  octave:4> A(1,1:2) = v'
  A =

    7  8  3
    4  8  8
\end{verbatim}


\subsection{Solving linear systems of equations}

\begin{itemize}
\item \verb+A\b+ solves the equation $Ax = b$.
\end{itemize}


\subsection{Inverses, decompositions, eigenvalues}

\begin{itemize}
\item \texttt{B = inv(A)} computes the inverse of $A$.
\item \texttt{[L,U,P] = lu(A)} computes the LU-decomposition
$LU=PA$.
\item \texttt{[Q,R] = qr(A)} computes the QR-decomposition $QR=A$.
\item \texttt{R = chol(A)} computes the Cholesky
decomposition of $A$.
\item \texttt{S = svd(A)} computes the singular values of $A$.
\item \texttt{H = hess(A)} brings $A$ into Hessenberg form.
\item \texttt{E = eig(A)} computes the eigenvalues of $A$.
\item \texttt{[V,D] = eig(A)} computes a diagonal matrix $D$,
containing the eigenvalues of $A$, and a matrix $V$ containing the
corresponding eigenvectors such that $AV=VD$.
\item \texttt{norm(X,p)} calculates the $p$-norm of vector $X$. If $X$
is a matrix, $p$ can only take the values \texttt{1}, \texttt{2} or
\texttt{inf}.  The default is $p=2$.
\item \texttt{cond(A)} computes the condition number of $A$ with
respect to the $2$-norm.
\end{itemize}
A lot of these commands support further options.  They can be listed
by typing \texttt{help funcname}.


\subsection{Testing for zero elements}

\begin{itemize}
\item \texttt{[i,j,v] = find(A)} finds the indices of all nonzero
elements of \texttt{A}.  The resulting vectors satisfy $A_{i(l),j(l)}
= v_l \neq 0$.
\item \texttt{any(v)} returns $1$ if the vector \texttt{v} contains
nonzero elements.
\item \texttt{any(A)} applies \texttt{any} to each of the columns of
the matrix \texttt{A}.
\end{itemize}


\section{Control structures}

\subsection{Functions}
\label{functions}

Traditionally, functions are also stored in plain text files with
suffix \texttt{.m}.  In contrast to scripts, functions can be called
with arguments, and all variables used within the function are
\emph{local}---they do not influence variables defined previously.

\paragraph*{Example:}

A function \texttt{f}, saved in the file named \texttt{f.m}.
\begin{verbatim}
  function y = f (x)
    y = cos(x/2)+x;
  end
\end{verbatim}

\paragraph*{Remark:}

In \Octave, several functions can be defined in a single script file.
\Matlab\, on the other hand, strictly enforces one function per
\texttt{.m} file, where the name of the function must match the name
of the file.  If compatibility with \Matlab\ is important, this
restriction should also be applied to programs written in \Octave.

\paragraph*{Example with two function values:}

A function \texttt{dolittle}, which is saved in the file named
\texttt{dolittle.m}.
\begin{verbatim}
  function [out1,out2] = dolittle (x)
    out1 = x^2;
    out2 = out1*x;
  end
\end{verbatim}
Calling the function:
\begin{verbatim}
  octave:1> [x1,x2]=dolittle(2)
  x1 = 4
  x2 = 8
  octave:2> whos
  *** currently compiled functions:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   user function         -      -  dolittle
  *** local user variables:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   real scalar           1      1  x1
   wd   real scalar           1      1  x2
\end{verbatim}

Obviously, the variables \texttt{out1} and \texttt{out2} were local to
\texttt{dolittle}.  Previously defined variables \texttt{out1} or
\texttt{out2} would not have been affected by calling
\texttt{dolittle}.


\subsection{Global variables}

\texttt{global name} declares \texttt{name} as a global variable.

\paragraph*{Example:}

A function \texttt{foo} in the file named \texttt{foo.m}:
\begin{verbatim}
  global N   % makes N a global variable; may be set in main file

  function out = foo(arg1,arg2)
  global N   % makes local N refer to the global N
  <Computation>
  end
\end{verbatim}
If you change \texttt{N} within the function, it changes in the
value of \texttt{N} everywhere.


\subsection{Loops}

The syntax of \emph{for}- and \emph{while}-loops is immediate from the
following examples:
\begin{verbatim}
  for n = 1:10
    [x(n),y(n)]=dolittle(n);
  end

  while t<T
    t = t+h;
  end
\end{verbatim}
\emph{For}-loop backward:
\begin{verbatim}
  for n = 10:-1:1
  ...
\end{verbatim}

\subsection{Branching}

Conditional branching works as follows.
\begin{verbatim}
  if x==0
    error('x is 0!');
  else
    y = 1/x;
  end

  switch pnorm
    case 1
      sum(abs(v))
    case inf
      max(abs(v))
    otherwise
      sqrt(v'*v)
  end
\end{verbatim}


\subsection{Functions of functions}

\begin{itemize}
\item \texttt{eval(string)} evaluates \texttt{string} as \Octave\ 
code.
\item \texttt{feval(funcname,arg1,$\dots$)} is equivalent to calling
\texttt{funcname} with arguments \texttt{arg1,$\dots$}.
\end{itemize}

\paragraph*{Example:}

Approximate an integral by the midpoint rule:
\[
  \int_a^b f(x) \, dx \approx \frac{b-a}{N} \sum_{j=0}^{N-1}
  f \left(a + (j+\frac{1}{2}) \frac{b-a}{N}\right)
\]
We define two functions, \texttt{gauss.m} and \texttt{mpr.m}, as
follows:
\begin{verbatim}
  function y = gauss(x)
    y = exp(-x.^2/2);
  end

  function S = mpr(fun,a,b,N)
    h = (b-a)/N;
    S = h*sum(feval(fun,[a+h/2:h:b]));
  end
\end{verbatim}
Now the function \texttt{gauss} can be integrated by calling:
\begin{verbatim}
octave:1> mpr('gauss',0,5,500)
\end{verbatim}


\subsection{Efficiency considerations}

Loops and function calls, especially through \texttt{feval}, have a
very high computational overhead. Therefore, if possible, vectorize
all operations.

\paragraph*{Example:}

We are programming the midpoint rule from the previous section
with a \emph{for}-loop (file name is \verb+mpr_long.m+):
\begin{verbatim}
  function S = mpr_long(fun,a,b,N)
    h = (b-a)/N; S = 0;
    for k = 0:(N-1),
      S = S + feval(fun,a+h*(k+1/2));
    end
    S = h*S;
  end
\end{verbatim}
We verify that \texttt{mpr} and \verb+mpr_long+ yield the same answer,
and compare the evaluation times.
\begin{verbatim}
  octave:1> t = cputime;
  > Int1=mpr('gauss',0,5,500); t1=cputime-t;
  octave:2> t = cputime;
  > Int2=mpr_long('gauss',0,5,500); t2=cputime-t;
  octave:3> Int1-Int2, t2/t1
  ans = 0
  ans = 45.250
\end{verbatim}


\subsection{Input and output}

\begin{itemize}
\item \texttt{save data var1 [var2 $\dots$]} saves the values of
variables \texttt{var1} etc.\ into the file \texttt{data}.
\item \texttt{load data} reads the file \texttt{data}, restoring the
values of the previously saved variables.
\item \texttt{fprintf(string[,var1,$\dots$])} resembles C syntax for
formatting output, see \texttt{man fprintf} under Unix.
\item \texttt{format [long|short]} enlarges or reduces the number of
decimals displayed.  Calling \texttt{format} without argument restores
the default.
\item \texttt{pause} Suspends evaluation until a key is pressed.
\end{itemize}

\paragraph*{Example:}
\begin{verbatim}
  octave:1> for k = .1:.2:.5,
  > fprintf('1/%g = %10.2e\n',k,1/k); end
  1/0.1 =   1.00e+01
  1/0.3 =   3.33e+00
  1/0.5 =   2.00e+00
\end{verbatim}


\section{Graphics}

\subsection{2D graphics}

\begin{itemize}
\item \texttt{plot(x,y[,fmt])} plots a line through the points
$(x_j,y_j)$. With the string \texttt{fmt} you can select line style
and color; see \texttt{help plot} for details.
\item \texttt{semilogx(x,y[,fmt])} like \texttt{plot} with a
logarithmic scale for the $x$ axis.
\item \texttt{semilogy(x,y[,fmt])} like \texttt{plot} with a
logarithmic scale for the $y$ axis.
\item \texttt{loglog(x,y[,fmt])} like \texttt{plot} with a logarithmic
scale for both axes.
\end{itemize}

Procedure for plotting a function $y=f(x)$:
\begin{enumerate}
\item Generate a vector with the $x$-coordinates of the points to be
plotted.
\begin{verbatim}
  x = x_min:step_size:x_max;
\end{verbatim}
(See also Section~\ref{vector}.)
\item Generate a vector containing the corresponding $y$-values by
letting $f$ act on the vector $x$ element-wise.
\begin{verbatim}
  y = f(x);
\end{verbatim}
\emph{Important:} Since $f$ operates element-wise, you must use of the
operators \texttt{.+}, \texttt{.-}, \texttt{.\^{}} etc.\ instead of
the usual \texttt{+}, \texttt{-} and \texttt{\^{}}!  (See
Section~\ref{element}.)
\item Finally you call the plot command.
\begin{verbatim}
  plot(x,y)
\end{verbatim}
\item To generate a coordinate grid, write:
\begin{verbatim}
  plot(x,y)
  grid
\end{verbatim}
\end{enumerate}


\paragraph*{Example:}

\begin{verbatim}
  octave:1> x = -10:.1:10;
  octave:2> y = sin(x).*exp(-abs(x));
  octave:3> plot(x,y)
  octave:4> grid
\end{verbatim}
\begin{center}
  \includegraphics[height=7cm]{2d-plot1} 
\end{center}

\subsection{3D graphics:}
\begin{itemize}
\item \texttt{[xx,yy] = meshgrid(x,y)} generates the grid data for a
3D plot---two matrices \texttt{xx} whose rows are copies of
\texttt{x}, and \texttt{yy} whose columns are copies of \texttt{y}.
\item \texttt{mesh(x,y,z)} plots a surface in 3D.
\end{itemize}

\paragraph*{Example:}
\begin{verbatim}
  octave:1> x = -2:0.1:2;
  octave:2> [xx,yy] = meshgrid(x,x);
  octave:3> z = sin(xx.^2-yy.^2);
  octave:4> mesh(x,x,z);
\end{verbatim}
\begin{center}
  \includegraphics[height=8cm]{3d-plot1} 
\end{center}


\subsection{Commands for 2D and 3D graphics}

\begin{itemize}
\item \texttt{title(string)} writes \texttt{string} as title for
the graphics.
\item \texttt{xlabel(string)} labels the $x$-axis with
\texttt{string}.
\item \texttt{ylabel(string)} labels the $y$-axis with
\texttt{string}.
\item \texttt{zlabel(string)} labels the $z$-axis with
\texttt{string}.
\item \texttt{axis(v)} set axes limits for the plot.  \texttt{v} is a
vector of form \texttt{v = (xmin, xmax, ymin, ymax[, zmin zmax])}.
\item \texttt{hold [on|off]} controls whether the next graphics output
should or not clear the previous graphics.
\item \texttt{clg} clears the graphics window.
% Only necessary in old version?!
%\item \texttt{replot} forces the current graphics window to be
%redrawn.  This is necessary in \Octave{} when altering an existing
%plots.  (This command exists only in \Octave{} and it is necessary
%only there.)
\end{itemize}


\section{Exercises}

\subsection{Linear algebra}

Take a matrix $A$ and a vector $b$ with
\[
A = \begin{pmatrix}
  1 & 2 \\
  3 & 4
\end{pmatrix} \quad \text{and} \quad
b = \begin{pmatrix}
  36 \\
  88
\end{pmatrix}.
\]
Solve the system of equations $Ax = b$.  Calculate the LU and QR
decompositions, and the eigenvalues and eigenvectors of $A$. Compute
the Cholesky decomposition of $A^T A$, and verify that
$\operatorname{cond}(A^T A) = \operatorname{cond}(A)^2$.

\paragraph*{Solution:}

\begin{verbatim}
  A = reshape(1:4,2,2).'; b = [36; 88];
  A\b
  [L,U,P] = lu(A)
  [Q,R] = qr(A)
  [V,D] = eig(A)
  A2 = A.'*A;
  R = chol(A2)
  cond(A)^2 - cond(A2)
\end{verbatim}

\subsection{Timing}

Compute the matrix-vector product of a $100 \times 100$ random matrix
with a random vector in two different ways.  First, use the built-in
matrix multiplication \texttt{*}.  Next, use \emph{for}-loops.
Compare the results and computing times.

\paragraph*{Solution:}

\begin{verbatim}
  A = rand(100); b = rand(100,1);
  t = cputime;
  v = A*b; t1 = cputime-t;
  w = zeros(100,1);
  t = cputime;
  for n = 1:100,
    for m = 1:100
      w(n) = w(n)+A(n,m)*b(m);
    end
  end
  t2 = cputime-t;

  norm(v-w), t2/t1
\end{verbatim}
Running this script yields the following output.
\begin{verbatim}
  ans = 0
  ans = 577.00
\end{verbatim}


\subsection{Stability functions of BDF-integrators}

Calculate all the roots of the polynomial
\[
  \frac{147}{60}\zeta^6 -6\zeta^5 + \frac{15}{2}\zeta^4
  - \frac{20}{3}\zeta^3 + \frac{15}{4}\zeta^2 - \frac{6}{5}\zeta +
  \frac{1}{6} .
\]
Hint: Use the command \texttt{compan}.

Plot these roots as points in the complex plane and draw a unit circle
for comparison. (Hint: \texttt{hold, real, imag}).

\paragraph*{Solution:}

\begin{verbatim}
  bdf6 = [147/60 -6 15/2 -20/3 15/4 -6/5 1/6];
  R = eig(compan(bdf6));
  plot(R,'+'); hold on
  plot(exp(pi*i*[0:.01:2]));
  if any(find(abs(R)>1))
    fprintf('BDF6 is unstable\n');
  else
    fprintf('BDF6 is stable\n');
  end
\end{verbatim}

\begin{center}
  \includegraphics[height=8cm]{bdf6} 
\end{center}


\subsection{3D plot}

Plot the graph of the function
\[
  f(x,y) = \exp(-x^2-y^2) .
\]

\paragraph*{Solution:}

\begin{verbatim}
  x = -3:0.1:3;
  [xx,yy] = meshgrid(x,x);
  z = exp(-xx.^2-yy.^2);
  figure, mesh(x,x,z);
  title('exp(-x^2-y^2)');
\end{verbatim}


\subsection{Hilbert matrix}

For each $n \times n$ Hilbert matrix $H$ where $n=1,\dots,15$ compute
the solution to the linear system $H x = b$, $b = $
\texttt{ones(n,1)}. Calculate the error and the condition number of
the matrix and plot both in semi-logarithmic coordinates.  (Hint:
\texttt{hilb, invhilb}.)

\paragraph*{Solution:}

\begin{verbatim}
  err = zeros(15,1); co = zeros(15,1);
  for k = 1:15
    H = hilb(k);
    b = ones(k,1);
    err(k) = norm(H\b-invhilb(k)*b);
    co(k) = cond(H);
  end
  semilogy(1:15,err,'r',1:15,co,'x');
\end{verbatim}


\subsection{Least square fit of a straight line}

Calculate the least square fit of a straight line to the points
$(x_j,y_j)$, given as two vectors $x$ and $y$.  Plot the points and
the line.

\paragraph*{Solution:}

\begin{verbatim}
  function coeff = least_square (x,y)
    n = length(x);
    A = [x ones(n,1)];
    coeff = A\y;
    plot(x,y,'x');
    hold on
    interv = [min(x) max(x)];
    plot(interv,coeff(1)*interv+coeff(2));
  end
\end{verbatim}


\subsection{Trapezoidal rule}

Write a program to integrate an arbitrary function $f$ in one variable
on an interval $[a,b]$ numerically using the trapezoidal rule with $h
= (b-a)/N$:
\[
  \int_a^b f(x) \, dx \approx h\left( \frac{f(a)+f(b)}{2} +
  \sum_{j=1}^{N-1} f(a+jh) \right).
\]
For a function $f$ of your choice, check, by generating a doubly
logarithmic error plot, that the trapezoidal rule is of order $2$.

\paragraph*{Solution:}

\begin{verbatim}
  function S = trapez(fun,a,b,N)
    h = (b-a)/N;
    % fy = feval(fun,[a:h:b]); better:
    fy = feval(fun,linspace(a,b,N+1));
    fy(1) = fy(1)/2;
    fy(N+1) = fy(N+1)/2;
    S = h*sum(fy);
  end

  function y = f(x)
    y = exp(x);
  end

  for k=1:15;
    err(k) = abs(exp(1)-1-trapez('f',0,1,2^k));
  end
  loglog(1./2.^[1:15],err);
  hold on;
  title('Trapezoidal rule, f(x) = exp(x)');
  xlabel('Increment');
  ylabel('Error');
  loglog(1./2.^[1:15],err,'x');
\end{verbatim}

\end{document}
