LV

1

Introduction to MATLAB
R. Gismondi
FIRM-Research Institute for Computational Methods
2
Index
3 10 51 64 86 99 111 118 127

• Appetizer of Numerical Computing
• Introduction
• Zeros and Roots
• Linear Equations, Norm and Condition Numbers
• Interpolation
  
• Quadrature
• Random Numbers
• Least Square
• Literature

Riccardo Gismondi
Folie 2
3

• Appetizer of Numerical Computing

Riccardo Gismondi
Folie 3
4
Appetizer of numerical computing (1)

• Let us start with a nice scientific (mathematical
  !!)
• computing experience
• It is possible to show ( one to everyone who
  brings the proof
• until next class ) that the sequence
  , with
• 
  converge to , i.e.
• .
• We are very good applied mathematicians, we
  have MATLAB,
• also we dont care about such easy task

Riccardo Gismondi
Folie 4
5
Appetizer of numerical computing (2)

• Let us write a MATLAB function (code) and let the
  computer
• work for us.(usually we are very lazy..). Let
  us try with n10

pitrue,estpi,relfel,absfelestimate_pi(10)
ans 3.141592653589793 2.000000000000000
36.338022763241867 1.141592653589793
3.141592653589793 2.828427124746190
9.968368384289397 0.313165528843603
3.141592653589793 3.061467458920717
2.550464159556757 0.080125194669076
3.141592653589793 3.121445152258050
0.641314885579486 0.020147501331743
3.141592653589793 3.136548490545931
0.160560696438413 0.005044163043862
3.141592653589793 3.140331156954739
0.040154685032539 0.001261496635054
3.141592653589793 3.141277250932756
0.010039578386341 0.000315402657037
3.141592653589793 3.141513801144145
0.002509951299956 0.000078852445648
3.141592653589793 3.141572940365566
0.000627491415994 0.000019713224227
3.141592653589793 3.141587725279960
0.000156872974192 0.000004928309833
Riccardo Gismondi
Folie 5
6
Appetizer of numerical computing (3)

• It looks very well, we try with n13

3.141592653589793 2.000000000000000
36.338022763241867 1.141592653589793
3.141592653589793 2.828427124746190
9.968368384289397 0.313165528843603
3.141592653589793 3.061467458920717
2.550464159556757 0.080125194669076
3.141592653589793 3.121445152258050
0.641314885579486 0.020147501331743
3.141592653589793 3.136548490545931
0.160560696438413 0.005044163043862
3.141592653589793 3.140331156954739
0.040154685032539 0.001261496635054
3.141592653589793 3.141277250932756
0.010039578386341 0.000315402657037
3.141592653589793 3.141513801144145
0.002509951299956 0.000078852445648
3.141592653589793 3.141572940365566
0.000627491415994 0.000019713224227
3.141592653589793 3.141587725279960
0.000156872974192 0.000004928309833
3.141592653589793 3.141591421504635
0.000039218488652 0.000001232085158
3.141592653589793 3.141592345611076
0.000009803267028 0.000000307978717
3.141592653589793 3.141592576545004
0.000002452411794 0.000000077044789
Riccardo Gismondi
Folie 6
7
Appetizer of numerical computing (4)

• We feel confident, the error (relative and
  absolute
• go to zero), we try with n30 ( bigger n, better
  results!!)

3.141592653589793 2.000000000000000
36.338022763241867 1.141592653589793
3.141592653589793 2.828427124746190
9.968368384289397 0.313165528843603
3.141592653589793 3.061467458920717
2.550464159556757 0.080125194669076
3.141592653589793 3.121445152258050
0.641314885579486 0.020147501331743
3.141592653589793 3.136548490545931
0.160560696438413 0.005044163043862
3.141592653589793 3.140331156954739
0.040154685032539 0.001261496635054
3.141592653589793 3.141277250932756
0.010039578386341 0.000315402657037
3.141592653589793 3.141513801144145
0.002509951299956 0.000078852445648
3.141592653589793 3.141572940365566
0.000627491415994 0.000019713224227
3.141592653589793 3.141587725279960
0.000156872974192 0.000004928309833
3.141592653589793 3.141591421504635
0.000039218488652 0.000001232085158
3.141592653589793 3.141592345611076
0.000009803267028 0.000000307978717
3.141592653589793 3.141592576545004
0.000002452411794 0.000000077044789
Riccardo Gismondi
Folie 7
8
Appetizer of numerical computing (5)
3.141592653589793 3.141592633463248
0.000000640647821 0.000000020126545
3.141592653589793 3.141592654807589
0.000000038763653 0.000000001217796
3.141592653589793 3.141592645321215
0.000000263197021 0.000000008268578
3.141592653589793 3.141592607375720
0.000001471039646 0.000000046214073
3.141592653589793 3.141592910939672
0.000008191701075 0.000000257349879
3.141592653589793 3.141594125195191
0.000046842654660 0.000001471605398
3.141592653589793 3.141596553704819
0.000124144517004 0.000003900115026
3.141592653589793 3.141674265021757
0.002597772561974 0.000081611431964
3.141592653589793 3.141829681889201
0.007544845100675 0.000237028299408
3.141592653589793 3.142451272494133
0.027330688571580 0.000858618904340
3.141592653589793 3.142451272494133
0.027330688571580 0.000858618904340
3.141592653589793 3.162277660168379
0.658424208974066 0.020685006578586
3.141592653589793 3.162277660168379
0.658424208974066 0.020685006578586
3.141592653589793 3.464101615137754
10.265779084358394 0.322508961547961
3.141592653589793 4.000000000000000
27.323954473516270 0.858407346410207
Riccardo Gismondi
Folie 8
9
Appetizer of numerical computing (6)

• Whats happened ???? Is the computer crazy ???
• Is the code false ???
• Or is Mr. Gismondi not such a good applied
  mathematician ???
• no job anymore by FIRM ???
• arbeitslos in Vienna ???
• Before to answer such bad questions, it is
  better, maybe,
• to learn what MATLAB is and what scientific
  computing means ()

() Such problems could have catastrophic
consequences, like 1. First Golf war
(1991) one Patriot-rocket exploded in a US
military point. 2. June, 4, 1996
ARIANE-rocket exploded
Riccardo Gismondi
Folie 9
10

• 1. Introduction

Riccardo Gismondi
Folie 10
11
What is MATLAB?

• A numerical computing environment
• A high-level programming language
• Matrix and vector-based
• Can be used for
• Mathematics and computation
• Easy vector and matrix entering and manipulation
• Plotting of information
• Implementation of algorithms
• Interfacing with programs in other languages

12
The command line

• MATLAB as a calculator
• gtgt (34)2 (ans 49)
• gtgt sin(90(3.14159/180)) (ans 1.0000)
• gtgt cos(pi) (ans -1)
• Assigning
• gtgt a4
• gtgt b2
• gtgt cab (c 16)

13
Variables

• General Variables
• ans (standard output)
• pi (ans 3.1416)
• i (ans 0 1.0000i)
• inf (ans Inf)
• New (Assigned) Variables
• a12
• b1sin(pi/2) (b1 1)
• csqrt(5a1-b1) (c 3)
• vvvwelcome (vvv welcome)

14
Operators (1)

• Arithmetical Operators
• (assignment)
• , - (addition, subtraction)
• (power function)
• /, \ (rigth, left division)
• (multiplication)

15
Operators (2)

• Logical Operators
• (logical not)
• (equal)
• lt, gt (less, greater)
• lt, gt (less or equal, greater or equal)
• , (logical and, logical or)

Riccardo Gismondi
Folie 15
16
Functions (1)

• General Functions
• who All current variables
• whos All current variables and their size
• clear a Clear variable (function) a
• clear all Clear all variables functions and
  debugging breakpoints
• clear Clear all variables and functions
• gwd The current directory
• cd Change the current directory
• format Set output format
• clc Clears all input and output from the Command
  Window

17
Functions (2)

• Help-Functions
• help elfun Help for all functions
• help elmat Help for all matrices
• help(f ) Help for the function f
• help ops Help for all Operators
• help plot Help for Plotting

18
Functions (3)

• Mathematic-Functions
• sin, cos, tan, cot
• exp, log, sqrt
• abs, angle, konj, imag, real
• fix, floor, ceil, round, sign
• fft, dft

19
Functions (4)

• Inline Functions
• Define a inline function
• finline('sin(x)/cos(x)')
• f
• f
• Inline function
• f(x) sin(x)/cos(x)
• f(pi/4) (ans 1.0000)

20
Vectors (1)

• Defining a Vector (1)
• Defining by square braces
• v1 7 3 4 5 6 (v 1 7 3 4 5
  6)
• v15 (v 1 2 3 4 5 )
• v10.52 (v 1.0000 1.5000 2.0000)
• Defining by assignment
• v1v(1)1v(2)22v(3)33v(4)43v(5)555
• v (v 1 22 33 43 555 )

21
Vectors (2)

• Defining a Vector (2)
• The transpose of v is defined by
• v
• ans
• 1
• 2
• 3
• 4
• 5

22
Vectors (3)

• Accessing Elements
• v15
• v(3) (ans 3)
• v(24) (ans 2 3 4)
• v(125) (ans 1 3 5)
• v(3)v(4)v(1) (ans 7)
• v(24)v(35) (ans 5 7 9)

Riccardo Gismondi
Folie 22
23
Vectors (4)

• Basic operations on vectors
• (assignment)
• , - (addition, subtraction)
• (power function)
• . (componentwise power function)
• (multiplication)
• . (componentwise multiplication)
• /, \ (rigth, left division)
• ./, .\ (componentwise right, left division)

Riccardo Gismondi
Folie 23
24
Vectors (5)

• Basic operations on vectors (Examples) (1)
• u15 (u 1 2 3 4 5)
• v5-11 (v 5 4 3 2 1)
• uv (ans 6 6 6 6 6)
• u-v (ans -4 -2 0 2 4)
• u.2 (ans 1 4 9 16 25)
• uv (ans 35)
• u.v (ans 5 8 9 8 5)
• u./v (ans 5.0000 2.0000 1.0000 0.5000
  0.2000)
• u.\v (ans 0.2000 0.5000 1.0000 2.0000
  5.0000)

Riccardo Gismondi
Folie 24
25
Vectors (6)

• Basic operations on vectors (Examples) (2)
• u/v (ans 0.6364)
• u\v
• ans
• 0 0 0 0 0
• 0 0 0 0 0
• 0 0 0 0 0
• 0 0 0 0 0
• 1.0000 0.8000 0.6000 0.4000 0.2000

Riccardo Gismondi
Folie 25
26
Matrices (1)

• Defining Matrices (1)
• Defining by square braces
• m 1 2 3 4 5 6 7 8 9
• m
• 1 2 3
• 4 5 6
• 7 8 9

Riccardo Gismondi
Folie 26
27
Matrices (2)

• Defining Matrices (2)
• Defining as a row of column vectors
• m 1 4 7' 2 5 8' 3 6 9'
• m
• 1 2 3
• 4 5 6
• 7 8 9

Riccardo Gismondi
Folie 27
28
Matrices (3)

• Defining Matrices (3)
• The transpose matrix of m is defined by
• m
• ans
• 1 4 7
• 2 5 8
• 3 6 9

Riccardo Gismondi
Folie 28
29
Matrices (4)

• Accessing Elements
• m 1 2 3 4 5 6 7 8 9
• m(2,3) (ans 6)
• m(12,23)
• ans
• 2 3
• 5 6
• m(2,2)m(1,1) (ans 6)

Riccardo Gismondi
Folie 29
30
Matrices (5)

• Basic operations on matrices
• (assignment)
• , - (addition, subtraction)
• (power function)
• . (componentwise power function)
• (multiplication)
• . (componentwise multiplication)
• /, \ (rigth, left division)
• ./, .\ (componentwise rigth, left division)

Riccardo Gismondi
Folie 30
31
Matrices (6)

• Basic operations on matrices (Examples) (1)
• m 1 2 3 4 5 6 7 8 9
• n 1 4 7 2 5 8 3 6 9
• mn
• ans
• 2 6 10
• 6 10 14
• 10 14 18

Riccardo Gismondi
Folie 31
32
Matrices (7)

• Basic operations on matrices (Examples) (2)
• mn
• ans
• 14 32 50
• 32 77 122
• 50 122 194
• n2
• ans
• 30 66 102
• 36 81 126
• 42 96 150

Riccardo Gismondi
Folie 32
33
Matrices (8)

• Basic operations on matrices (Examples) (3)
• n.m
• ans
• 1 8 21
• 8 25 48
• 21 48 81
• m./n
• ans
• 1.0000 0.5000 0.4286
• 2.0000 1.0000 0.7500
• 2.3333 1.3333 1.0000

Riccardo Gismondi
Folie 33
34
Programming in MATLAB (1)

• If Statement
• Syntax
• if expression statements end

Riccardo Gismondi
Folie 34
35
Programming in MATLAB (2)

• If-else Statement(1)
• Syntax
• if expression statements1else
  statements2end

Riccardo Gismondi
Folie 35
36
Programming in MATLAB (3)

• switch Statement (1)
• Syntax
• switch switch_expr
• case case_expr
• statement, ..., statement
• case case_expr1, case_expr2, case_expr3,
  ...
• statement, ..., statement
• otherwise
• statement, ..., statement
• end

Riccardo Gismondi
Folie 36
37
Programming in MATLAB (4)

• switch Statement (2)
• Example
• method 'Bilinear'
• switch lower(method)
• case 'linear','bilinear'
• disp('Method is linear')
• otherwise
• disp('Unknown method.')
• end
• Method is linear

Riccardo Gismondi
Folie 37
38
Programming in MATLAB (5)

• while loop
• Syntax
• while expression statements end

a1b4 while altb aa1 end a 2 a 3 a
4
Riccardo Gismondi
Folie 38
39
Programming in MATLAB (6)

• for loop
• Syntax
• for xinitvalendval statements end
• for xinitvalstepvalendval statements end

Riccardo Gismondi
Folie 39
40
Programming in MATLAB (8)

• M-File Scripts
• Have no input or output arguments
• Variable are stored in MATLAB workspace
• Usually implemented in MATLAB workspace
• M-File Functions
• Are program routines
• Accept input arguments and return output
  arguments
• Usually implemented in M-files

Riccardo Gismondi
Folie 40
41
Programming in MATLAB (9)

• Script (Example)
• Fibonacci Numbers
• M-file function to compute the first 10
  Fibonacci Numbers
• Definiton f(1)1, f(2)1, f(n)f(n-1)f(n-2)
• fzeros(1,10) f(1)1 f(2)1
• for ii310
• f(ii)f(ii-1)f(ii-2)
• end
• f (f 1 1 2 3 5 8 13 21
  34 55)

Riccardo Gismondi
Folie 41
42
Programming in MATLAB (10)

• Declare M-file function (1)
• Syntax
• function out1, out2, ... funname(in1,in2,
  ...)
• Online Help
• Comments, eg, Input, Output,
• Function code
• Defines function funname with inputs in1, in2,
  etc. and returns outputs out1, out2, etc.

Riccardo Gismondi
Folie 42
43
Programming in MATLAB (11)

• Declare M-file function (2)
• Example Fibonacci Numbers
• function f fibonacci(n) An M-file function to
  compute the Fibonacci Numbers
• Definiton f(1)1, f(2)1,
  f(n)f(n-1)f(n-2)
• fzeros(n,1)f(1)1 f(2)2
• for kk3n f(kk)f(kk-1)f(kk-2)
• end
• f(10) (f 1 1 2 3 5 8 13
  21 34 55 )
• 55

Riccardo Gismondi
Folie 43
44
Graphics and Plotting (1)

• Basic Plotting
• plot plot a a function or a data
• title outputs the title.
• xlabel, ylabel labels the x-axis, y-axis of the
  current axes
• grid grid lines for 2-D and 3-D plots
• legend outputs a legend.
• hold retain current graph in figure
• axis manipulates commonly used axes properties

Riccardo Gismondi
Folie 44
45
Graphics and Plotting (2)

• Basic Plotting (Example) (1)
• x -pi0.1pi
• plot(x,sin(x)) title('sin(x)')xlabel('x'),
  ylabel('y')

Riccardo Gismondi
Folie 45
46
Graphics and Plotting (3)

• Basic Plotting (Example) (2)
• X1,X2meshgrid(-pi.05pi)
• Z((pi2)/4)(sin(piX1).sin(piX2))
• figure plot3(X1,X2,Z)
• legend('f(x_1,x_2)\pi2/4sin(\pix_1)sin(\pix_2)
  ')
• xlabel('x_1','fontsize',12,'fontweight','b')
• ylabel('x_2','fontsize',12,'fontweight','b')
• zlabel('f(x_1,x_2)','fontsize',12,'fontweight','b'
  )

Riccardo Gismondi
Folie 46
47
Graphics and Plotting (4)

• Basic Plotting (Example) (3)

48
Graphics and Plotting (5)

• Mesh and Surface
• meshgrid generate X and Y arrays for 3-D plots
• mesh draws a wireframe mesh
• surf plot a 3-D shaded surface
• colorbar add a colorbar tool
• patch creating patch graphics objects
• hidden remove hidden lines from mesh plot

Riccardo Gismondi
Folie 48
49
Graphics and Plotting (7)

• Mesh and Surface (Example) (1)
• X1,X2meshgrid(-pi.05pi)
• Z((pi2)/4)(sin(piX1).sin(piX2))
• figure surf(X1,X2,Z)
• legend('f(x_1,x_2)\pi2/4sin(\pix_1)sin(\pix_2)
  ')
• xlabel('x_1','fontsize',12,'fontweight','b')
• ylabel('x_2','fontsize',12,'fontweight','b')
• zlabel('f(x_1,x_2)','fontsize',12,'fontweight','b'
  )

Riccardo Gismondi
Folie 49
50
Graphics and Plotting (8)

• Mesh and Surface (Example) (2)

51

• 2. Zeros and Roots

Riccardo Gismondi
Folie 51
52
Bisection

• Computing of
• Try x 1,5. This x is too big
• Try x 1,25. This x is too small
• Try x 1.375. This x is to small
• Approximations to
• 1,5 1,25 1,375 1.3125 1.4063 ..

Riccardo Gismondi
Folie 52
53
Bisection

• Computing of
• MATLAB function
• function app root2()M2a1b2eps0.001k0
  while b-a gt eps app(ab)/2 if
  app2gtM bapp else
  aapp end kk1display(app)end

Riccardo Gismondi
Folie 53
54
Bisection

• Solving f(x)0 with Bisection
• Find a very small interval a, b, on which f
  changes sign.
• Divides the interval in two by computing
  c(ab)/2
• There are now two possibilities
• f(a) and f(c) have opposite signs
• f(c) and f(b) have opposite signs
• Applied bisection to the sub-interval where the
  sign change occurs

Riccardo Gismondi
Folie 54
55
Bisection

• Solving f(x)0 with Bisection
• MATLAB function
• function app bisection(f,a,b,eps)kk0while
  b-a gt epsabs(b) app(ab)/2
  if sign(f(app)) sign(f(b)) bapp
  else aapp
  endkkkk1 display(app)end

Riccardo Gismondi
Folie 55
56
Bisection

• Solving f(x)0 with Bisection
• MATLAB function (testing)
• finline('x2-2')
• app bisection(f,0,2,10-10)
• app(end) 1.414213562267833
• sqrt(2) 1.414213562373095

Riccardo Gismondi
Folie 56
57
Newtons Method

• Solving f(x)0 with Newton Method
• Newtons Iterations
• with a starting value

Riccardo Gismondi
Folie 57
58
Newtons Method

• Solving f(x)0 with Newton Method
• The MATLAB function
• function app newton(f,x,eps,h,met)xprev0
  qq0while abs(x-xprev)gt epsabs(x)
  xprevx if met 1
  fder(f(xh)-f(x))/h elseif met2
  fder (f(xh)-f(x-h))/(2h) end
• xx-f(x)/fder appx qqqq1
  display(app)
• end

Riccardo Gismondi
Folie 58
59
Newtons Method

• Solving f(x)0 with Newton Method
• MATLAB function (testing)
• finline('x2-2')
• app newton(f,0.01,10-5,10-5,1)
• app(end) 1.414213562378965
• sqrt(2) 1.414213562373095

Riccardo Gismondi
Folie 59
60
Newtons Method

• Solving f(x)0 with Newton Method
• MATLAB function (testing a perverse example)
• finline('sign(x-2)sqrt(abs(x-2))')
• newton(f,0.01,10-5,10-5,1)
• app
• .

Riccardo Gismondi
Folie 60
61
Secant Method

• Solving f(x)0 with Secant Method
• Iterations
• with a starting values

Riccardo Gismondi
Folie 61
62
Secant Method

• Solving f(x)0 with Secant Method
• The MATLAB function
• function app secant(f,a,b,eps)k0while
  abs(a-b)gt epsabs(b) ca ab
  bb(b-c)/(f(c)/f(b)-1) appb
  kk1display(app)end

Riccardo Gismondi
Folie 62
63
Secant Method

• Solving f(x)0 with Secant Method
• MATLAB function (testing)
• finline('x2-2')
• secant(f,0.01,0.02)
• app
• 66.6733 0.0500 0.0799 15.4271 0.2085
  0.3336 3.8177 0.7886 1.0878 1.5231 1.4006
  1.4137 1.4142
• ans
• 1.4142

Riccardo Gismondi
Folie 63
64

• 3. Linear Equations, Norm and
• Condition Numbers

Riccardo Gismondi
Folie 64
65
Solving Linear Systems

• A quadratic system of linear equations has form
• Ax bwhere A is a given (n, n)-Matrix and b is
  a given n-Vector
• The solution can be written
• x A b
• Example

Riccardo Gismondi
Folie 65
66
The MATLAB Backslash Operator

• System of equations (1)
• AX B
• A matrix of any size, B matrix with as many rows
  as A
• Solution (MATLAB)
• X A\B
• System of equations (2)
• XA B
• A matrix of any size, B matrix with as columns
  rows as A
• Solution (MATLAB)
• X B/A

Riccardo Gismondi
Folie 66
67
The MATLAB Backslash Operator

• System of equations (3)
• Example (MATLAB) (1)
• A8 1 6 3 5 7 4 9 2
• B1 2 3
• XA\B
• X
• 0.0500
• 0.3000
• 0.0500

Riccardo Gismondi
Folie 67
68
The MATLAB Backslash Operator

• System of equations (4)
• Example (MATLAB) (2)
• A1 2 3 4 5 6 7 8 9
• B1 2 3
• XA/B
• X
• 1.0000
• 2.2857
• 3.5714

Riccardo Gismondi
Folie 68
69
Permutation and Triangular Matrices

• Permutation matrix
• Identity matrix with rows and columns
  interchanged.
• Exactly one 1 in each row and column
• All other elements are 0
• Example
• P
• Multiplying a matrix A on the left by permutation
  matrix permutes the rows of A.
• Multiplying on the right permutes the columns of
  A

Riccardo Gismondi
Folie 69
70
Permutation and Triangular Matrices

• Permutation matrix
• Example (MATLAB)
• A1 2 3 4 5 6 7 8 9
• p3 1 2
• A(p,)
• ans
• 7 8 9
• 1 2 3
• 4 5 6

Riccardo Gismondi
Folie 70
71
Permutation and Triangular Matrices

• Triangular Matrices
• Upper triangular matrix
• Has all its nonzero elements above or on the main
  diagonal
• Unit lower triangular matrix
• Has ones in the main diagonal
• All the rest of its nonzero elements below the
  main diagonal

Riccardo Gismondi
Folie 71
72
Permutation and Triangular Matrices

• Triangular Matrices
• Example
• U
  L
• Linear equations involving triangular matrixes
  are also trivial to solve

Riccardo Gismondi
Folie 72
73
LU Factorization

• LU factorization for a Matrix A
• LU PAP permutation matrix, L unit lower
  triangular, U upper triangular.
• exists for all nonsingular matrix A and is
  unique.
• The system Ax b become
• LUx Pb
• Solve lower triangular system
• Ly Pb
• Finally, solve upper triangular system
• Ux y

Riccardo Gismondi
Folie 73
74
LU Factorization

• Example (MATLAB)
• A 1 2 3 4 5 6 7 8 1
• b6 15 16
• The solution to Ax b is obtained with matrix
  division
• xA\b
• x
• 1
• 1
• 1

Riccardo Gismondi
Folie 74
75
LU Factorization

• Example (MATLAB)
• Solution to Ax b by LU-Factorization
• L,U,P lu(A)
• L
• 1.0000 0 0
• 0.1429 1.0000 0
• 0.5714 0.5000 1.0000

Riccardo Gismondi
Folie 75
76
LU Factorization

• Example (MATLAB)
• Solution to Ax b by LU-Factorization
• U
• 7.0000 8.0000 1.0000
• 0 0.8571 2.8571
• 0 0 4.0000
• P
• 0 0 1
• 1 0 0
• 0 1 0

Riccardo Gismondi
Folie 76
77
LU Factorization

• Example (MATLAB)
• Solution to Ax b by LU-Factorization
• Verify that LU is a permuted version of A
• PA LU
• ans
• 1 1 1
• 1 1 1
• 1 1 1

Riccardo Gismondi
Folie 77
78
LU Factorization

• Example (MATLAB)
• Solution to Ax b by LU-Factorization
• y L\Pb
• y
• 16.0000
• 3.7143
• 4.0000

Riccardo Gismondi
Folie 78
79
LU Factorization

• Example (MATLAB)
• Solution to Ax b by LU-Factorization
• x U\y
• x
• 1.0000
• 1.0000
• 1.0000

Riccardo Gismondi
Folie 79
80
Norm and Condition Numbers

• Norm of a Vector x
• p1 Manhattan norm
• p2 Euclidean norm
• Chebyshev norm
  

Riccardo Gismondi
Folie 80
81
Norm and Condition Numbers

• Norm of a Vector (MATLAB)
• norm(x,p) computes
• norm(x) computes
• norm(x,inf) computes

Riccardo Gismondi
Folie 81
82
Norm and Condition Numbers

• Norm of a Vector (MATLAB)
• Example
• x14norm1norm(x,1)norm2norm(x)norm3norm(x,i
  nf)
• Produces
• x 1 2 3 4
• norm1 10
• norm2 5.4772
• norm3 4

Riccardo Gismondi
Folie 82
83
Norm and Condition Numbers

• Condition Number of Matrix A
• Condition Number of Matrix (MATLAB)
• cond(A,1) computes
• cond(A) computes
• cond(A,inf) computes

Riccardo Gismondi
Folie 83
84
Norm and Condition Numbers

• Condition Number of Matrix (MATLAB)
• Example
• A2 3 7 1 10 4 3 3 3
• cond(A)ans 7.0453
• cond(A,inf)ans 10.4762
• cond(A,1)ans 9.6508

Riccardo Gismondi
Folie 84
85
Norm and Condition Numbers

• Condition Number of Matrix (MATLAB)
• Example (bad condition)
• A1 1/2 1/3 1/2 1/3 1/4 1/3 1/4 1/5
  (Hilbert Matrix)
• cond(A)ans 524.0568
• cond(A,inf)ans 748.0000
• cond(A,1)ans 748.0000

Riccardo Gismondi
Folie 85
86

• 4. Interpolation

Riccardo Gismondi
Folie 86
87
The Interpolating Polynomial

• What is interpolation?
• Given a set of n data points n1
• with
• for
• There is a unique polynomial p of degree less
  than n with
• called nodes
• is called a data point
• p is called interpolant

Riccardo Gismondi
Folie 87
88
The Interpolating Polynomial

• Lagrange interpolation
• with

Riccardo Gismondi
Folie 88
89
The Interpolating Polynomial

• Lagrange interpolation (MATLAB function
  la_interp)
• v la_interp(x,y,u) finds L(x) with L(x(j))
  y(j) and returns v(k) L(u(k)).
• function v la_interp(x,y,u)nlength(x)vzeros
  (size(u))for k 1n wones(size(u))
  for j 1k-1 k1n w(u-x(j))./(x(k)-x
  (j)).w end vvwy(k)end

Riccardo Gismondi
Folie 89
90
The Interpolating Polynomial

• Lagrange interpolation (testing function
  la_interp) (1)
• x03
• y-5 -6 -1 16
• u-.25 .1 3.25
• vla_interp(x,y,u)

Riccardo Gismondi
Folie 90
91
The Interpolating Polynomial

• Lagrange interpolation (testing function
  la_interp) (2)
• v Columns 1 through 13 -4.5156 -4.7034
  -4.9001 -5.0999 -5.2966 -5.4844 -5.6571
  -5.8089 -5.9336 -6.0254 -6.0781 -6.0859
  -6.0426
• Columns 14 through 26-5.9424 -5.7791
  -5.5469 -5.2396 -4.8514 -4.3761 -3.8079
  -3.1406 -2.3684 -1.4851 -0.4849 0.6384
  1.8906
• Columns 27 through 363.2779 4.8061
  6.4814 8.3096 10.2969 12.4491 14.7724
  17.2726 19.9559 22.8281

Riccardo Gismondi
Folie 91
92
The Interpolating Polynomial

• Lagrange interpolation (testing function
  la_interp) (2)
• plot(x,y,'o',u,v,'-')

Riccardo Gismondi
Folie 92
93
Piecewise Linear Interpolation

• What is Piecewise Linear Interpolation?
• Connect straight lines between data points
• Any intermediate value read off from straight
  line
• Interpolant
• Linear function that passes through
  and

Riccardo Gismondi
Folie 93
94
Piecewise Linear Interpolation

• MATLAB function piecelin
• function v piecelin(x,y,u)PIECELIN Piecewise
  linear interpolation. v piecelin(x,y,u) finds
  piecewise linear L(x) with L(x(j)) y(j) and
  returns v(k) L(u(k)).
• First divided difference
• delta diff(y)./diff(x)

Riccardo Gismondi
Folie 94
95
Piecewise Linear Interpolation

• MATLAB function piecelin
• Find subinterval indices k so that x(k) lt u
  ltx(k1)n length(x)k ones(size(u))for j
  2n-1 k(x(j) lt u) jend
• Evaluate interpolants u - x(k)v y(k)
  s.delta(k)

Riccardo Gismondi
Folie 95
96
Piecewise Linear Interpolation

• MATLAB testing function piecelin (1)
• x16
• y16 18 21 17 15 12
• u-.25 .1 3.25
• vpiecelin(x,y,u)

Riccardo Gismondi
Folie 96
97
Piecewise Linear Interpolation

• MATLAB testing function piecelin (2)
• v Columns 1 through 1313.5000 13.7000
  13.9000 14.1000 14.3000 14.5000 14.7000
  14.9000 15.1000 15.3000 15.5000 15.7000
  15.9000
• Columns 14 through 2616.1000 16.3000 16.5000
  16.7000 16.9000 17.1000 17.3000 17.5000
  17.7000 17.9000 18.1500 18.4500 18.7500
• Columns 27 through 3619.0500 19.3500
  19.6500 19.9500 20.2500 20.5500 20.8500
  20.8000 20.4000 20.0000

98
Piecewise Linear Interpolation

• MATLAB testing function piecelin (3)
• plot(x,y,'o',u,v,'-')

Riccardo Gismondi
Folie 98
99

• 5. Quadrature

Riccardo Gismondi
Folie 99
100
Adaptive Quadrature

• Let be
• We seek to compute
• If m is a point between a and b than

Riccardo Gismondi
Folie 100
101
Adaptive Quadrature

• Idea
• The value of the integral is the area under curve

Riccardo Gismondi
Folie 101
102
Adaptive Quadrature

• Idea
• If we can approximate each of the two integrals
  within a specified tolerance, then the sum gives
  us the result.
• If not, we can recursively apply the additive
  property to each of the intervals a, xi and
  xi, b.

Riccardo Gismondi
Folie 102
103
Basic Quadrature Rules

• Midpoint rule
• Approximation
• with
  and

Riccardo Gismondi
Folie 103
104
Basic Quadrature Rules

• Trapezoid rule
• Approximation
• with

Riccardo Gismondi
Folie 104
105
Basic Quadrature Rules

• Midpoint and Trapezoid rule (Example)

Riccardo Gismondi
Folie 105
106
Basic Quadrature Rules

• Midpoint and Trapezoid rule (Example)
• We get the errors
• S M 1/3 - 1/4 1/12
• S T 1/3 - 1/2 -1/6
• If the error in T were -2 times the error in M
  than
• S T - 2(S M)
• S 2M/3 T/3
• S is the exact value of the integral
• In any case, S is usually a more accurate
  approximation than either M or T alone.

Riccardo Gismondi
Folie 106
107
Basic Quadrature Rules

• Simpsons Rule

Riccardo Gismondi
Folie 107
108
Basic Quadrature Rules

• MATLAB Example(1)
• Midpoint, Trapezoid and Simpson's approximation
  forthe function sqrt(1-x2) for two subintervals
• a 0 b 1 h(b-a) m 1/2
• y0 sqrt(1-a2) y1 sqrt(1-m2) y2
  sqrt(1-b2)
• M 2hy1 mid-point ruleT h(y0y2)/2
  trapezoidal ruleS h(y04y1y2)/3
  Simpson ruleIquad('sqrt(1-x.2)',0,1) MATLAB

Riccardo Gismondi
Folie 108
109
Basic Quadrature Rules

• MATLAB Example(2)
• Midpoint Trapezoid and Simpson's approximation
  forthe function sqrt(1-x2) for two subintervals
• M 1.732051
• T 0.5000
• S 1.488034
• I 0.7854

Riccardo Gismondi
Folie 109
110
Numerical Integration (MATLAB)

• Commands
• quad(fun,a,b)
• quad(fun,a,b,tol)
• Examples
• quad('x.2',0,1) ans 0.3333
• quad('sqrt(1-x.2) ',0,1,10-8)
• ans 0.7854

Riccardo Gismondi
Folie 110
111

• 6. Random Numbers

Riccardo Gismondi
Folie 111
112
Pseudorandom numbers

• Random numbers (MATLAB)
• By start up a fresh MATLAB
• format long randans 0.814723686393179
• The first MATLAB random number
• All MATLAB users keep getting this same number
• True random? No.
• Computers are (in principle) deterministic
  machines and should not exhibit random behavior.

Riccardo Gismondi
Folie 112
113
Pseudorandom numbers

• What is a random sequence of numbers?
• by professor D. H. Lehmer
• A random sequence is a vague notion in which
  each term is unpredictable to the uninitiated and
  whose digits pass a certain number of tests
  traditional with statisticians

Riccardo Gismondi
Folie 113
114
Pseudorandom numbers

• MATLAB
• Y rand (returns a pseudorandom)
• Y rand(n) returns an (n, n)-matrix of values
  derived as described above.
• Y rand(m,n) (returns an random (m, n)-matrix)

Riccardo Gismondi
Folie 114
115
Pseudorandom numbers

• MATLAB (Example)
• YY 0.8147
• Y rand(3)Y 0.9058 0.6324 0.5469
  0.1270 0.0975 0.9575 0.9134 0.2785
  0.9649Y
• Y rand(2,3)Y 0.4218 0.7922 0.6557
  0.9157 0.9595 0.0357

Riccardo Gismondi
Folie 115
116
Uniform Distribution

• Algorithm to generate random numbers
• a, c and m initial values ( integers)
• MATLAB function
• function u_rand rand r_gen(a, c, m, x)
  random numbers generatorxvecxfor j2m5
  m5 random integers x rem((axc),m)
  (ax c) mod m xvec(j)x randxvec
  u_randrand/m
• end

Riccardo Gismondi
Folie 116
117
Uniform Distribution

• Algorithm to generate random numbers
• MATLAB function (testing)
• r_gen(13,0,31,1)
• ans Columns 1 through 2113 14 27 10
  6 16 22 7 29 5 3 8
  11 19 30 18 17 4 21 25
  15Columns 22 through 369 24 2 26
  28 23 20 12 1 13 14 27
  10 6 16

Riccardo Gismondi
Folie 117
118

• 7. Least square

Riccardo Gismondi
Folie 118
119
Models and Curve Fitting

• Curve Fitting
• Given a set of m data points
• with
• for
• We seek to find a function y such that

Riccardo Gismondi
Folie 119
120
Models and Curve Fitting

• Curve Fitting (Solution)
• Idea
• model y by a linear combination of n basis
  functions
• In Vector-Matrix notation
• The least squares solution (MATLAB)

Riccardo Gismondi
Folie 120
121
Models and Curve Fitting

• Curve Fitting (Models)
• Linear
• Polynomial

Riccardo Gismondi
Folie 121
122
Models and Curve Fitting

• Curve Fitting (MATLAB)
• polyfit - computes least squares polynomial fits
• Example
• xd 00.55 yd sin(xd) grad 5p1
  polyfit(xd,yd,grad) p1p1
• -0.0051 0.0823 -0.3728 0.2277 0.9119
  0.0019

Riccardo Gismondi
Folie 122
123
Models and Curve Fitting

• Curve Fitting (MATLAB)
• x 00.15 y1 polyval(p1,x)plot(xd,yd,'ro',x
  ,y1,'b-')

Riccardo Gismondi
Folie 123
124
Householder Reflections

• What is Householder Reflection?
• Matrix of the form
• with
• H reflects every vector x in the hyper plane
  span(u) .
• H is both symmetric and orthogonal
• H is never formed
• The application of H to a vector x

Riccardo Gismondi
Folie 124
125
The QR Factorization

• What is QR Factorization?
• Given a (n, m)-Matrix A, with m gt n.
• QR Factorization for A
• Q is (m, m)-Matrix and orthogonal, and R is (n,
  n)-Matrix and upper triangular
• Can be used to solve linear systems (see LU
  Factorization), least squares problems, etc.

Riccardo Gismondi
Folie 125
126
The QR Factorization

• MATLAB Example
• A1 2 3 4 5 6 7 8 9
• Q,Rqr(A)
• Q -0.1231 0.9045 0.4082 -0.4924
  0.3015 -0.8165 -0.8616 -0.3015 0.4082
• R -8.1240 -9.6011 -11.0782 0
  0.9045 1.8091 0 0
  -0.0000

Riccardo Gismondi
Folie 126
127

• 8. Literature

Riccardo Gismondi
Folie 127
128
Literature

• 1. Numerical Computing with Matlab
• http//www.mathworks.com/moler/chapters.html
  
• 2. Wissenschftliches Rechnen mit MATLAB
• Quarteroni, Saleri, Springer, 2005
• 3. Applied Numerical Methods for Engineers
  Using
• MATLAB and C
•   Schilling, Harris, Thomson Engineering,
  2000
• 4. Introduction to MATLAB 7 for Engineers
• Palm, McGraw-Hill, 2005

Riccardo Gismondi
Folie 128