Transcript MA/CS 375 - Rice U - Computational and Applied Mathematics

MA/CS 375
Fall 2002
Lecture
Summary Week 1 Week 7

Theoretical Exercise
• Say we wish to design a self-playing computer
game like asteroids.
• The player controls a rocket.
• There are enemy rockets who can shoot
torpedoes and who can ram the player.
• There are moving asteroids of fairly arbitrary
shape.
• Let’s work through the details involved……
Review
• We are going to do a fast review from
basics to conclusion.
Week 1
• Basics
Ok – Cast your minds back
We are first faced with a text prompt:
Basic Math
We can create our own variables and apply basic
algebraic operations on these variables.
Basic operations include addition (+), multiplication (*),
subtraction (-), division (/ or \).
a=1
b = a/2
c = b*3
Do the math and:
b = a/2 = ½
c = b*3 = 1.5
Basic Arrays (Matrices)
• We can create one- or two-dimensional
variables of arbitrary size.
• We can then apply basic linear algebra
operations on these.
Example 2
A is a matrix with
3 rows and 2 columns.
 2.1 3.23 


A  4.12 1.893


 7.1

5 

Matrix Addition in Matlab
2
A
3
4
B
4
1
2 
2
1 
CAB
Matrix Subtraction in Matlab
2
A
3
4
B
4
1
2 
2
1 
C AB
Matrix Multiplication
• There is a specific definition of matrix
multiplication.
C  A B
• In index notation:
Cij 
NcolsA
A
ik
Bkj
k 1
• i.e. for the (i,j) of the result matrix C we take the
i’th row of A and multiply it, entry wise, with the
j’th column of B
Example 4
(matrix multiplication)
2
A
3
4
B
4
1

2
2
1 
C  A B
Functions in Matlab
• Matlab has a number of built-in functions:
– cos, sin, tan, acos, asin, atan
– cosh, sinh, tanh, acosh, asinh, atanh
– exp, log, sqrt
• They all take matrices as arguments and
return matrices of the same dimensions.
• e.g. cos([1 2 3])
• For other functions type: > help matlab\elfun
Example of Function of Vector
Custom-Made Matlab
Functions
Say we wish to create a function that turns Cartesian coordinates into
polar coordinates. We can create a text file with the following text. It
can be called like any built in function.
function [ radius, theta] = myfunc( x, y)
% this is a comment, just like // in C++
% now create and evaluate theta (in radians)
theta = atan2(y,x);
% now create and evaluate radius
radius = sqrt( x.^2 + y.^2);
function end
Custom Built Function For Matlab
• Make sure you use cd
in Matlab to change to the
directory containing
myfunc.m
• the arguments to myfunc
could also have been
matrices – leading to two
matrices being output.
Matlab as Programming Language
• We can actually treat Matlab as a coding
language.
• It allows script and/or functions.
• Loops are allowed, but since Matlab is an
interpreted language, their use can lead to
slow code.
Loops in Matlab
One variant of Matlab loop syntax is:
for var=start:end
commands;
end
Example of a Matlab Loop
• Say I want to add the
numbers from 1 to
10, without using the
Matlab intrinsic sum.
Week 2
• Plotting
• Finite precision effects
Plotting
• Recall we can create a function of a
vector.
• Then plot the vector as one ordinate and
the function of the vector as the other
ordinate.
Example 5
(figure created by Matlab)
Adding Titles, Captions, Labels,
Multiple Plots
subplot
• If you wish to create a figure with two
sub-figures you can use the subplot
function:
• subplot(1,2,1) requests
1 row of figures
2 columns of figures
1st figure
subplot(1,2,1)
subplot(1,2,2)
Starting Numerics
• We next considered some limitations
inherent in fixed, finite-precision
representation of floating point numbers.
A Convergent Binary
Representation of Any Number
Between 0 and 1
n 
1
m   mn n where mn  1 or 0
2
n 1
a similar representation in base 10:
n 
1
d   d n n where d n  0,1,2,3,4,5,6,7,8 or 9
10
n 1
Volunteer ?
Finite Binary Approximations of a
Real Number
n N
1
m   mn n +TN where mn  1 or 0
2
n 1
1
We can easily show that TN is bounded as: TN  N
2
(think of Zeno’s paradox)
Monster #1
• Consider:
• What should its behavior be as:
f  x
1 x 1


x 0
x
• Use subplots for the following 2 plots
• Plot this function at 1000 points in:
x   1,1
• Plot this function at 1000 points in:
 4 ,4 
• Label everything nicely, include your name in the title.
• In a text box explain what is going on, print it out and hand it
in
f  x
1 x 1


x 0
Monster #1
x
Each stripe is a region
where 1+ x is a constant
(think about the gaps between
numbers in finite precision)
Then we divide by x and the
stripes look like hyperbola.
The formula looks like
(c-1)/x with a new c for each
stripe.
when we zoom in we see that the large+small operation is introducing
order eps errors which we then divide with eps to get O(1) errors !.
Monster #2
• Consider:
f  x   e log(1  e )
x
x
• What should its behavior be as: x  
Monster #2 cont
(finite precision effects from large*(1+small) )
As x increases past ~=36 we
see that f drops to 0 !!
Limit of

f  x   e x log(1  e  x )
lim f  x   lim e log 1  e
x 
x 
x
x

 log 1  e  x  
 lim 


x
x  

e


 e  x 
 1  e x 
 lim 
by l'Hopital's rule


x
x 
  e 


 1 
 lim 
1


x
x  1  e


Monster #4
• Consider: g  x   cos  x    sin  x     sin  x  





• What should its behavior be as:   0

Monster 4 cont
Behavior as delta  0 :
 sin  x     sin  x  
 sin  x  cos    cos( x )sin    sin  x  
lim 
 lim 


 0
 0






 sin  x   cos    1  cos( x )sin   
 lim 

 0 



 cos  x 
or if you are feeling lazy use the definition of derivative, and remember:
d(sin(x))/dx = cos(x)
Monster 4 cont
(parameter differentiation, delta=1e-4)
OK
Monster 4 cont
(parameter differentiation, delta=1e-7)
OK
Monster 4 cont
(parameter differentiation, delta=1e-12)
Worse
Monster 4 cont
(parameter differentiation, delta=1e-15)
Bad
When we make the delta around about machine precision we see
O(1) errors !.
Approximate Explanation of
Monster #4
f  x     f  x    f '  x   O  2 
1) Taylor’s thm:
yˆ  xˆ  ˆ
 x    O ( )
2) Round off errors
3) Round off in computation
of f and x+delta
fˆ  yˆ   f  yˆ   O ( )
 f  x     O ( )
fˆ  yˆ   fˆ  xˆ   f  y   f  x   O  
4) Put this together:
  f '  x   O  2   O   
ˆf  yˆ   fˆ  xˆ    f '  x   O  2   O  
 fˆ  yˆ   fˆ  xˆ  
 

 O    O  









i.e. for 

or equivalently 

  10
8
approximation error decreases as delta
decreases in size.
BUT for 
  10
8
round off dominates!.
Week 3 & 4
• We covered taking approximate
derivatives in Matlab and manipulating
images as matrices.
Week 5
• Approximation of the solution to ordinary
differential equations.
• Adams-Bashforth schemes.
• Runge-Kutta time integrators.
Ordinary Differential Equation
• Example:
u 0  a
du
 u
dt
u T   ?
• t is a variable for time
• u is a function dependent on t
• given u at t = 0
• given that for all t the slope of
us is –u
• what is the value of u at t=T
Forward Euler
Numerical Scheme
• Numerical scheme:
u
n
 u

t
Discrete scheme:
u
•
n 1
u
n 1
n
 1  t  u
n
where: un  approximate solution at t=nt
Summary of dt Stability
• 0 < dt <1
stable and convergent since as
dt  0 the solution approached
the actual solution.
• 1 <= dt < 2 bounded but not cool.
• 2 <= dt exponentially growing,
unstable and definitely not cool.
Application: Newtonian Motion
N-Body Newtonian
Gravitation Simulation
• Goal: to find out where all the objects are after
a time T
• We need to specify the initial velocity and
positions of the objects.
• Next we need a numerical scheme to advance
the equations in time.
• Can use forward Euler…. as a first approach.
Numerical Scheme
For m=1 to FinalTime/dt
For n=1 to number of objects
v
m 1
n
 v  dt
m
n
x
iN

i 1,i  n
m
i
-x
x
End
End
 x  dt  v
m
n
m
n

M G
x x
m
i
End
For n=1 to number of objects
m 1
n
m
n
i
m 3
n 2
AB Schemes
Linear:
y n1  y n   t  f n
Quadratic:
y n1  y n 
Cubic:
3t n t n1
f 
f
2
2
Essentially we use
interpolation and a
Newton-Cotes
quadrature formula to
formulate:
23t n 16t n1 5t n2
y y 
f 
f 
f
12
12
12
Quartic:
55t n 59t n1 37 t n2 9t n3
n 1
n
y y 
f 
f 
f 
f
24
24
24
24
Quintic:
1901t n 2774t n1 2616t n2 1274t n3 251t n4
y n1  y n 
f 
f 
f 
f 
f
720
720
720
720
720
n 1
n
Runge-Kutta Schemes
• See van Loan for derivation of RK2 and
RK4.
• I prefer the following (simple) scheme
due to Jameson, Schmidt and Turkel
(1981):
y  yn
for k  s : 1:1
t
n
y=y 
f  y
k
end
y n 1  y
Runge-Kutta Schemes
• Beware, it only works when f is a function
of y and not t
y  yn
for k  s : 1:1
t
y=y 
f  y
k
n
end
y
n 1
y
• here s is the order of the scheme.
Week 6
• Colliding disks project
Week 7
• Matrix inverse and what can go wrong.
• Solving lower and upper triangular
systems.
• LU factorization and partial pivoting.
Recall
• Given a matrix M then if it exists the
inverse of a matrix M-1 is that matrix which
satisfies:
1 0
0 1
M 1M  MM 1  


0 0
0
0



1
Examples
 2 4
• If A  

1 3
• If
1
3
A
5

2
what is A 1
2 3 4
1 3 2

2 1 4

4 1 3
what is A 1
• If A is an NxN matrix how can we
calculate its inverse ? 
Example cont
• The inverse of A can be calculated as:
1 
  1
2





1
A 

1 
1 
 1

 2 
3
 
 
• Now let’s see how well this exact solution
works in Matlab:
1  
A 1

 
 
Example cont
1  
A 1

 
 

1 

1
  1
2



1
A 
1
 1
 2
3

 





With this formulation of
1
A the product of A and
its inverse only satisfies
the definition to 6 decimal
places for delta=0.001