Ch7.3 Part 2: Using Octave

To illustrate how to use basic matrix-vector commands in Octave, we will consider the sytem \( A\mathbf{x} = \mathbf{b} \) below.

\[ \begin{array}{rrrr} x_1 & + & 2x_2 & = & 5 \\ 3x_1 & - & 4 x_2 & = & 6 \end{array} \, \, \Leftrightarrow \begin{bmatrix} 1 & 2 \\ 3 & -4 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 5 \\ 6 \end{bmatrix} \]

In augmented matrix form,

\[ \begin{array}{rrrr} x_1 & + & 2x_2 & = & 5 \\ 3x_1 & - & 4 x_2 & = & 6 \end{array} \, \, \Leftrightarrow \begin{bmatrix} 1 & 2 & 5 \\ 3 & -4 & 6 \end{bmatrix} \]

Who invented the round table?

Sir Cumference

The following code shows how to enter the augmented matrix \( [A | \mathbf{b}] \) in Octave.

Ab = [3,2,5;3,-4,6]

Ab =

   3   2   5
   3  -4   6

\[ \begin{array}{rrrr} x_1 & + & 2x_2 & = & 5 \\ 3x_1 & - & 4 x_2 & = & 6 \end{array} \, \, \Leftrightarrow \begin{bmatrix} 1 & 2 & 5 \\ 3 & -4 & 6 \end{bmatrix} \]

To find the reduced row echelon form of the augmented matrix, use the rref command:

Ab = [3,2,5;3,-4,6];
rref(Ab)

ans =

   1.0000        0   1.7778
        0   1.0000  -0.1667

The LU decomposition of \( A \) can be found as follows:

A = [3,2;3,-4];
[L,U]=lu(A)

L =

   1   0
   1   1

U =

   3   2
   0  -6

• The following code shows how to enter the augmented matrix in Python.
• As you can see, this is more difficult in Python than in Octave, and would be problematic for large matrices.

import numpy as np
Ab = np.array([[3,2,5],[3,-4,6]])
print(Ab)

[[ 3  2  5]
 [ 3 -4  6]]

• Python does not have an rref command.
  
• The LU decomposition of \( A \) can be found as follows:

from scipy.linalg import lu
A = np.array([[3,2],[3,-4]])
lu(A)

(array([[1., 0.],
       [0., 1.]]), array([[1., 0.],
       [1., 1.]]), array([[ 3.,  2.],
       [ 0., -6.]]))

• Many of our numerical methods use for loops.
• These are also a little easier in Octave than in Python.
• The following code chunk shows how the range of the index k is interpreted in an Octave for loop.

for k = 1:4
k = k
end

k = 1
k = 2
k = 3
k = 4

• Indentation not required, but it helps visual organization.
• Note that the for loop finishes with the end command.

for k = 1:4
  k = k
end

k = 1
k = 2
k = 3
k = 4

If we don't want to output the values to the screen, then a semicolon at the end of the line will suppress output:

for k = 1:4
  k = k;
end

Without semicolons:

for k = 1:2
  k = k
end

k = 1
k = 2

• The lambda function that Python uses is replaced by the anonymous function in Octave, indicated by the @ symbol.
• Both functions and vectors use round parentheses instead of parentheses and square brackets in Python.
• The Octave program is written entirely as a script, as opposed to definitions in Python, because markdown implementation in Octave doesn't allow for Octave functions.

f = @(x) exp(-x^2);
for k = 1:n  
  nodes(k) = a + (k-1)*h;   
  yvalues(k) = f(nodes(k)); 
end  

#Initial conditions
n = 10;   #semicolon at end to suppress output
a = 0;    #Left endpoint of [a,b]
b = 1;    #right endpoint of [a,b] 
h = (b-a)/n;   #Compute step size h  

#Initialize vectors. 
#zeros(M,1) returns Mx1 matrix of zeros. 
nodes = zeros(n+1,1);  
yvalues = zeros(n+1,1);

#Function to integrate
f = @(x) exp(-x^2);

#For n rectangles, there are n+1 nodes. 
for k = 1:n  #Includes 1 and n
  nodes(k) = a + (k-1)*h;   #left nodes
  yvalues(k) = f(nodes(k)); #left yvalues
end     

#Initialize S and add up yvalues
S = 0;
for k = 1:n
S = S + yvalues(k);
end

#Left sum value
S = h*S

#Compute "actual" value and PRE. 
Q = quad(f,0,1)
PRE = (Q-S)/Q*100

The output when running script is the same as for Python.

S = 0.7778
Q = 0.7468
PRE = -4.1499

• As for the left sum method, the Octave program is written as a script rather than as a function.
• In Octave markdown, plots are printed to a window.
  
  • Plot doesn't appear in output unless we print it to a file and then display the file.

#Euler's Method
for n = 1:N
  x(n+1) = x(n)+h;              #Next x value
  y(n+1) = y(n)+h*f(x(n),y(n)); #Next y value
  gn(n+1) = g(x(n+1));          #Next g value
  sn(n+1) = n+1;                #Next s value
end
#Print plot to file using this filename
print("7.3Fig1.png") 

#Initial conditions
N = 20;
a = 0;    #Left endpoint of [a,b]
b = 1;    #right endpoint of [a,b] 
h = (b-a)/N;   #Compute step size h  
x0 = a;        #Initial x value as given in IVP
y0 = 0;        #Initial y value as given in IVP

#Derivative function
f = @(x,y) x + y;

#Analytic (exact) solution of IVP for this problem. 
g = @(x) exp(x) - x - 1;

#Initialize vectors
x = zeros(N+1,1);  #node vector
y = zeros(N+1,1);  #numerical solution vector 
gn = zeros(N+1,1); #exact solution vector
sn = zeros(N+1,1); #index (subscript) vector

#Initial conditions
x(1) = x0;
y(1) = y0;
gn(1) = g(x0);
sn(1) = 1;    

#Euler's Method
for n = 1:N
  x(n+1) = x(n)+h;              #Next x value
  y(n+1) = y(n)+h*f(x(n),y(n)); #Next y value
  gn(n+1) = g(x(n+1));          #Next g value
  sn(n+1) = n+1;                #Next s value
end

#Display results as columns of matrix R. 
#Use space key to position headers
#The "%5.4f" specifies decimal format ("f"), with
#5 digits total, 4 to right of decimal point.

R = [sn x y gn]';
fprintf('n  x(n)  y(n)  Exact \n');  
fprintf('%2d %2.2f %5.4f  %5.4f \n', R); 

#Plot numerical solution
#Set up new figure and clear previous figure

figure(1); clf(1) 
plot(x(1:N+1),y(1:N+1),'bo')
hold all #Retain for multiple plots

L = N*h;  #[x0,x0+L] = x-axis interval
xvalues = linspace(x0,x0+L,1000);  

plot(xvalues,g(xvalues),'r') #Plot exact solution

#Plot style commands
axis([0 1.1 0 0.8]) #Viewing window; no commas
xlabel('x')
ylabel('y')
k = legend('Euler Method','Analytic Solution', 'Location','northwest');
set(k,'fontsize',12);
set(gca,'fontsize',12)  #gca = "get current axis" being used

#Print plot to file using this filename
print("7.3Fig1.png")  

Here are the numerical results for \( N=10 \).

n  x(n)  y(n)  Exact 
 1 0.00 0.0000  0.0000 
 2 0.10 0.0000  0.0052 
 3 0.20 0.0100  0.0214 
 4 0.30 0.0310  0.0499 
 5 0.40 0.0641  0.0918 
 6 0.50 0.1105  0.1487 
 7 0.60 0.1716  0.2221 
 8 0.70 0.2487  0.3138 
 9 0.80 0.3436  0.4255 
10 0.90 0.4579  0.5596 
11 1.00 0.5937  0.7183 

The graph shown below is for \( N=20 \).

• Octave is a free clone of MATLAB.
• MATLAB is widely used in STEM fields.
• MATLAB stands for “Matrix Laboratory” and is designed for working with matrix algebra and systems of linear equations.
• A student version of MATLAB costs between $49 and $100 depending on the packages purchased.
  

• Octave is easier to use than Python for matrix-vector applications.
• The programming style used for our numerical methods is very similar to what we did in Python, but there are small differences in the commands and notation used.
• We will use Octave within a markdown environment.
• Greater versatility is available directly within Octave.