How to clean a python script from annoying ^M ?
Lec 13. Roots of Equations (Bracketing Methods)
The problem of finding the real roots (zeros) of a given continuous function arises frequently
in science and engineering.
• Roots of a quadratic equation ax2 + bx + c = 0 are given by
√
−b ± b2 − 4ac
x1,2 =
2a
• How about the cubic equation ax3 +bx2 +cx+d = 0? Or, nonlinear equations involving
different types of functions!
A. Graphical Approach
A simple method for obtaining the estimate of the root of the equation f (x) = 0 is to make
a plot of the function and observe where it crosses the x-axis.
• Consider the cubic equation x3 + x2 − 3x − 3 = 0.
a. First, of course, we ask: Are there any? Plot it to get an idea what the curve
y = x3 + x2 − 3x − 3 looks like. Try Gnuplot.
$ gnuplot
gnuplot> plot x**3+x**2-3*x-3
It crosses the x-axis and looks very flat on it, …many roots? (think: at most how
many roots could there be for this f ? ) Also plot y = 0 (x-axis) for better display.
gnuplot> plot 0, x**3+x**2-3*x-3
Notice the very big range on the y-axis which makes small variations invisible, so
we can’t really see much at this scale.
b. To resolve it better, let’s plot it on a smaller interval, say, [-2,2]:
gnuplot> plot [-2:2]
0, x**3+x**2-3*x-3
It looks much better, doesn’t it? It crosses the axis three times, so there are three
real roots. Let’s magnify the scale (zoom in) around the largest root (inside [1,2]),
so try:
gnuplot> plot [1:2]
0, x**3+x**2-3*x-3
Now it’s better resolved, and we can zoom in further:
gnuplot> plot [1.6:1.8]
0, x**3+x**2-3*x-3
c. Repeat the zooming in until you can estimate the root to three decimals, and write
it down. This is the Graphical Method for finding roots of (simple) functions! Can
you “guess” the exact value of the root ? Write it down.
B. Bisection (Interval Halving) Method
We can also estimate an interval for the root using the Intermediate Value Theorem for
continuous functions (function changes sign on opposite sides of a root). That is,
If f (x) = 0 and a ≤ x ≤ b, then f (a)f (b) < 0.
To find a root xr of the equation f (x) = 0 accurate within a specified tolerance , we choose
an interval [xl , xu ] such that f (xl )f (xu ) < 0 and use the following algorithm.
Repeat
Set xr = 21 (xl + xu )
If f (xr ) = 0 then xr is the zero.
Stop
If f (xl )f (xr ) < 0 then
Set xu = xr
Else if f (xl )f (xr ) > 0 then
Set xl = xr
End if
Until (|xu − xl | < 2)
Absolute Error. Let L0 = ∆x0 = x0u − x0l be the length of the initial interval (step 0).
Then the absolute error in each iteration can be computed as
Ea0 = L0
L0
2
L0
Ea2 = 2
2
..
.
Ea1 =
L0
2n
• The final value of xr , obtained after n bisections, approximates the root and it is in
error by not more than L0 /2n . The method is slow, but it is guaranteed for convergence.
Ean =
• The number of bisections n needed to approximate the root xr within a desired error
tolerance Ea,d can be estimated “a priori” by
L0
ln(L0 /Ea,d )
n = log2
=
Ea,d
ln 2
Examples.
1. Do three iterations of bisection method by hand to find a zero of f (x) = x6 − 2 on [0, 2].
Give an upper bound for the error.
(0)
Iteration 1: xl
(0)
(0)
(0)
(0)
(0)
(0)
(1)
(1)
(1)
(1)
(2)
(2)
(2)
(2)
=
2−0
= 0.25
23
= 0, xu = 2; f (xl )f (xu ) < 0; xr = 12 (xl + xu ) = 1.
Iteration 2: f (xl )f (xr ) > 0; xl
(0)
(0)
(1)
= 1, xu = 2; xr = 12 (xl + xu ) = 23 .
(1)
(1)
(2)
= 1, xu = 32 ; xr = 12 (xl + xu ) = 45 .
Iteration 3: f (xl )f (xr ) < 0; xl
Upper bound for the error:
(0)
Ea3 =
(0)
∆x0
xu − xl
=
3
2
23
2. If the desired absolute magnitude of the error Ea,d = 10−4 and the length of the interval
L0 = ∆x0 = 2, how many bisections will you have to do to get the required accuracy in the
solution?
n=
ln(2 × 104 )
ln(L0 /Ea,d )
=
= 14.3 ≈ 15
ln 2
ln 2
Lec 15. Modules in Python
Python has a large library of modules and packages that extends its
functionality.
I The import statement makes reference to the modules that are
Python files containing definitions and statements.
I If you write import in a program, the Python interpreter
executes the statements in the file .py and enters the
module name into the current namespace, so that the
functions/procedures it defines are available with the “dotted
syntax”: .. Recall how you used math
module:
import math
use math.sqrt(2), math.pi
from math import *
What does this do?
I You can also define your own module by placing a code within a file
module name.py somewhere in your computer (for small projects,
usually in the same directory where you run your code).
Note: DO NOT name your module anything that isn’t a valid
Python identifier (hyphen, starting with a digit, name of a built-in
functions, modules, etc.)
Lec 15. Numerical Integration
n
X
Z b
f (x)dx = lim
n→∞
a
f (xi∗ )∆x
i=1
I Rectangle Rule: (left sum, right sum, middle sum)
Z b
f (x)dx ≈
a
n
X
f (xi−1 )∆x,
i=1
I Trapezoid Rule:
n
X
f (xi )∆x,
i=1
Z b
a
n
X
(x̄i )∆x
i=1
"
#
n−1
X
∆x
f (x0 ) + 2
f (xi ) + f (xn )
f (x)dx ≈
2
i=1
xi−1 +xi
∆x = b−a
n , xi = a + i∆x, x̄i =
2
Next
Classwork 6.
Modules in Python: Numerical Integration
Date: 07/24/2022
2
√
π
Z x
0
2
e −t dt = erf(x)
Lec 16. Introduction to NumPy
NumPy is the fundamental package for scientific computing with
Python and is the base for many other packages.
I NumPy has a powerful custom N-dimensional array object for
efficient and convenient representation of data.
I NumPy supports vectorization: a single operation can be carried out
on an entire array, no need to loop over the array’s element!
I The absence of explicit looping and indexing makes the code
cleaner, less error-prone and closer to the standard mathematical
notation it reflects.
I NumPy has tools for integration with other programming languages
used for scientific programming like C/C++ and Fortran.
I Moreover, it will introduce you to a whole other set of libraries such
as SciPy and Scikit-learn, which you can use to solve almost any
problem in scientific computing and data science.
NumPy: Basic array methods
The most powerful construct of NumPy is the ndarray.
I The NumPy array class ndarray provides a generic container for
multidimensional array of a homogeneous data which can be sorted,
reshaped, written to and read from files and much more.
I That is, unlike Python lists and tuple, the element in NumPy array
cannot be of different types. The type in the ndarray is specified
by (optional) data type object (dtype).
I Import NumPy package with the statement import numpy as np
and then refer to its attributes with the prefix np. Try:
import numpy as np
x=np.array( [1,2,3] )
y=np.array( [[1, 2, 3],[4,5,6]]] )
z=np.array([[[1,2],[3,4]],[[5,6],[7,8]]])
print(x)
print(y)
print(z)
print(x[0],x[2], x[-1])
print(y[0,1],y[1,2],z[0,0,0],z[1,1,1])
print(type(y))
print(y.shape)
print(y.ndim)
print(y.dtype)
print(y.size)
NumPy: Basic array methods
I x.ndim returns the dimension of the array x.
I x.shape returns a tuple giving the number of elements in each
dimension of the array x. For example, the output (3,) means an
array of 3 elements, (2,3) means 2 rows with 3 columns, (2,2,2)
means 2 rows of arrays of 2 rows and 2 columns.
I x.size returns the total number of elements in the array x.
I x.dtype returns the data type of the elements in the array x. The
data type is “upcast” to the most general type if they are mixed but
compatible types. You can also set the data type in the np.array
construct. Try:
import numpy as np
x=np.array( [1,2.,3] )
print(x)
print(x.dtype)
x=np.array( [1, 2, 3+4j] )
print(x)
print(x.dtype)
x=np.array( [1,2,3],dtype=complex)
print(x.dtype)
print(x.nbytes)
# float, float32, float16, int, int8, int16, int32
NumPy: Creating/Initializing an array
I np.arange(n) creates a list of 0 to n-1 (like range, but it can
generate floating point sequence!)
I np.linspace(a,b,n) creates a sequence (1D array) of n evenly
spaced numbers between a and b. np.linspace has a couple of
optional boolean arguments. Setting retstep=True, returns the
step-size. endpoint=False omits the final point in the sequence.
Try:
import numpy as np
x=np.arange(10) ; y=np.arange(10,dtype=float)
print(x); print(y)
x=np.linspace(-np.pi,np.pi,20) #20 evenly spaced numbers between -pi and pi
x,dx=np.linspace(-np.pi,np.pi,20,retstep=True) #optional retstep=True returns the step-size
print(x); print(dx)
y=np.linspace(0,np.pi,10,endpoint=False)
print(y)
I Also try: np.geomspace(a,b,n) and np.logspace(a,b,n)
import numpy as np
u=np.geomspace(1,64,5)
print(u)
v=np.logspace(3,4,5)
print(v)
NumPy: Creating/Initializing an array
I You can create an array of a desired shape using np.empty. The
initial element values are undefined.
I np.zeros and np.ones create an array of specified shape with
elements pre filled with 0 and 1 respectively.
I np.empty, np.zeros and np.ones also take optional dtype
argument. The default data type is float.
I You can create new arrays of the same shape of the exiting arrays,
use np.empty like, np.zeros like or np.ones like. Try:
import numpy as np
x=np.empty((2,3)) # not initialized
y=np.zeros((3,2)) # prefilled with 0’s
z=np.ones((3,3))
# prefilled with 1’s
print(x); print(y); print(z)
x=np.ones((2,3), dtype=int)
print(x)
a=np.ones_like(y)
print(a)
b=np.arange(6).reshape(2,3)
print(b)
NumPy: Array operations
NumPy provides many familiar mathematical functions that the math
module does, implemented as universal functions that act on each
element of the array producing an array in return.
I Thus, no need to loop explicitly over the array elements. Try these
operations in one dimensional arrays, print each statement:
import numpy as np
x=np.linspace(-5,5,4)
xp1=x+1
xt2=x*2
xm1=x-1
x2=x**2
sqrtx=np.sqrt(x)
sinx=np.sin(x)
ex=np.exp(x)
y=np.zeros(4)+2
xmy=x-y
xoy=x/y
xty=x*y
I Also try the following:
import numpy as np
M=np.identity(4)
# creates an identity matrix of size 4x4
N=M+3
P=2*M
Q=M-N
R=M-x
NumPy: Reading and writing array to a file
In scientific computing, we need to read data frequently from a text file.
NumPy provides several functions for reading data from a text file.
I np.save saves NumPy arrays in (plateform-independent) binary file
format (.npy) that can be read (on any other operating system)
using the function np.load. Try:
M=np.random.rand(100,3)
np.save(’array1.npy’,M)
N=np.load(’array1.npy’)
print(M); print(N)
# saves the array M in a file array1.npy
# loads the array from array1.npy as N
I np.savetxt saves the data in human readable text-files.
np.loadtxt reads the data from text-files. Try:
np.savetxt(’array2.txt’,M)
P=np.loadtxt(’array2.txt’)
print(P)
I np.loadtxt has several optional arguments to handle the data file
that contains descriptive header and columns are aligned in a fixed
width format or separated by one or more delimiting characters
(spaces, tabs or commas). Syntax:
np.loadtxt( fname, dtype=, comments=’#’, delimiter=None, skiprows=0,
usecols=None, unpack=False)
NumPy: Reading and writing array to a file
The description of the arguments of np.loadtxt
I fname: only required argument, name of the file to read data from
I dtype: data type of the array, default is float
I comments: comments in a file usually (# or %) to be ignored
I delimiter: string used to separate columns in a file, default is
None, meaning any amount of whitespace (spaces, tabs) delimits
the data. To read comma-separated (csv) file, set delimiter=’,’
I skiprows: an integer giving the number of lines at the start of the
file to skip over before reading the data, default is 0 (no header)
I usecols: a sequence of column indexes determining which columns
of the file to return as data, default is None (all columns)
I unpack: by default data is returned as a single array of rows and
columns, unpack=True will transpose this array so that individual
columns can be picked off an assigned to different variables.
Revisit Classwork 5 after NumPy
In CW5 you learned how to read data from a file using f.readline(s)
and sys.stdin.readline(s). Here we use np.loadtxt for the same.
I First, go to the working directory MA305/CLASSWORK/CW4/. Make a
copy of cw5a.py as cw5np.py. Comment out all the lines related to
reading from the file.
I Import NumPy package with the statement import numpy as np.
I The file ’cw5a.txt’ has data in four columns, and you want 1st, 3rd
and 4th columns as D, X and Y. Also, you do not need to save the
first line (header)! Call the NumPy function np.loadtxt as:
D, X, Y = np.loadtxt(’cw5.txt’, skiprows=1, usecols=(0,2,3), unpack=True)
print()
I Print D, X and Y to check if it read the data. What did you
observe?
I There are two complains: n is not defined, D is not integer. Do the
following:
D, X, Y = np.loadtxt(’cw4a.txt’, skiprows=1, usecols=(0,2,3), unpack=True)
n=len(D)
D=D.astype(int) #this will change the data type of D float (default) to integer
Dr. H. Khanal
MA 305.01, Summer B, 2022
Final Project
Presentation: Thursday, 08/11/2022
t
Final Report: 11:59pm Saturday, 08/13/2022
! Choose one problem and work in the assigned group. Organize your python code using modules and
functions. Write sufficient comments in your code for clarity. Use formatted output for nicer display. Give
appropriate title, labels and legends in the plots. Prepare your final report using the guidelines posted at your
course Canvas, inserting tables and figures in the text with appropriate captions.
I
I. Computing ⇡.
1. Numerical Integration. The value of ⇡ can be computed by using the definite integrals.
Z 1
Z 1
4
1
p
(i)
dx
(ii)
dx
2
1
+
x
1
x2
0
1
(a) Implement the Trapezoid rule:
"
#
Z b
N
X1
1
1
f (x) dx ⇡ x f (x0 ) +
f (xi ) + f (xN ) ,
2
2
a
i=1
x=
b
a
N
and the Simpson’s 1/3rd rule:
2
3
Z b
N
N
X1
X2
x4
f (x)dx ⇡
f (x0 ) + 4
f (xi ) + 2
f (xi ) + f (xN )5 ,
3
a
i=1,3,5
i=2,4,6
,
x=
xi = a + i x
b
a
N
, xi = a + i x
to approximate the integral (i).
(b) Implement the Midpoint rule:
Z b
a
f (x)dx ⇡
x
N
X
f (xi ),
x=
b
i=1
a
N
, xi = a +
2i
1
2
x
to approximate both the integrals (i) and (ii).
Which method gives the best approximation to ⇡? How accurate did you get? Display your results
(error of approximation for different values of N ) both graphically and in tabular form.
2. Sum of Alternating Series. From Calculus II we know that
tan 1 x
=
x
x3
x5
+
3
5
This is an alternating series which converges for
···
⇡
N
X
n=1
( 1)n+1
x2n 1
2n 1
(1)
1 < x 1. Thus, for x > 0 we know that if we only
sum the first N terms then the error (the difference between arctan x and the partial sum) is no larger
than the absolute value of the next term, i.e |x2N +1 /(2N + 1)|.
1
a. Write a function arctan(x,N) that takes x and N , and returns the sum of the first N terms in
the series. The value of ⇡ can be approximated by 4*arctan(1,N). Explore how large should N
be to get a good approximation to ⇡!
b. Machin’s Formula
16 tan 1 ( 15 )
1
4 tan 1 ( 239
)=⇡
uses a combination of the series of arctan(x) to compute ⇡ . Write a function machine(N) that
1
takes N and returns the value 16*arctan( 15 ,N) – 4*arctan( 239
,N) as an approximation to ⇡.
How big should N be to get a good approximation to ⇡ using the function machine()? Any
surprise?
c. Compare your results from (a) and (b) with the Madhava’s series
✓
◆
p
1
1
1
12 1
+
+ ··· = ⇡
3 · 3 5 · 32
7 · 33
3. Dart Board Algorithm. The value of ⇡ can also be calculated using the following procedure (Dart
Board Algorithm): Inscribe a unit circle x2 + y 2 = 1 in a square (|x| 1, |y| 1). Randomly generate
N points in the square. Determine the number of points in the square that are also in the circle. Let
r be the number of points in the circle divided by the number of points in the square. Then ⇡ ⇡ 4r.
Write a function dartboard(N) to implement the Dart Board Algorithm. How many correct digits did
you manage to find?
2
II. Root Finding using Newton’s Method.
1. Write a function newton() to implement Newton’s method for finding a root of a given function F (x).
x0 = given;
xn+1 = xn
F (xn )
F 0 (xn )
n = 0, 1, 2, · · · , Nmax
A good way to decide convergence is to test both the relative change and the residual:
|1
xn+1 /xn | < T OL,
|F (xn )| < T OL
The main program should read in initial guess for the root x0 , and tolerance for testing convergence ✏ (these
could also be set in the code, if you prefer). Your program should output the input values (x0 and ✏) at the
beginning, then at each iteration print out: n, xn , F (xn ) and, upon convergence, the number of iterations
n, the root xn and the residual |F (xn )|. The values of F (xn ) and F 0 (xn ) may be computed in a function or
two separate functions.
a. Set T OL = 10 12 and find the real root of the “Newton’s cubic” x3
2x
5 = 0 (the only equation
Newton ever bothered to solve with his method, in 1671 !)
b. Use it to find the three roots of the cubic (from Lab 4) x3 + x2
3x
3 = 0. Is Newton method more
efficient than Bisection?
c. Commissioner Gordon has been found dead in his office. At 8:00pm, the county coroner determined
the core temperature of the corpse to be 90 F. One hour later, the core temperature had dropped
to 85 F. Maintenance reported that the building’s air conditioning unit broke down at 4:00pm. The
2
temperature in the commissioner’s office was 68 F at that time. The computerized climate control
system recorded that the office temperature rose at a rate of 1 F per hour after the air conditioning
stopped working. Captain Furillo believes that the infamous Doc B killed the commissioner. Doc B,
however, claims that he has an alibi. Lois Lane was interviewing him at the Daily Planet Building,
just across the street from the commissioner’s office. The receptionist at the Daily Planet Building
checked Doc B into the building at 6:35pm, and the interview tapes confirm that Doc B was occupied
from 6:40pm to 7:15pm. Could Doc B have killed the commissioner?
Hint: Assume that the core temperature of the corpse was 98.6 F at the time of death and began decreasing
immediately obeying Newton’s Law of Cooling. Let T (t) denote the temperature ( F) of the corpse at time t (hr).
Set t = 0 to correspond to 8:00pm. Then
dT
=
dt
k(T
72
t), T (0) = 90, T (1) = 85
Now, to answer the question you need to determine the exact time of death. That is, find t such that T (t) = 98.6 F.
3
III. Newton Method, Fractals and Julia Sets
1. Write a function newton() to implement Newton’s method for finding a root of a given function F (x).
x0 = given;
xn+1 = xn
F (xn )
F 0 (xn )
n = 0, 1, 2, · · · , Nmax
A good way to decide convergence is to test both the relative change and the residual:
|1
xn+1 /xn | < T OL,
|F (xn )| < T OL
The main program should read in initial guess for the root x0 , and tolerance for testing convergence ✏ (these
could also be set in the code, if you prefer). Your program should output the input values (x0 and ✏) at the
beginning, then at each iteration print out: n, xn , F (xn ) and, upon convergence, the number of iterations
n, the root xn and the residual |F (xn )|. The values of F (xn ) and F 0 (xn ) may be computed in a function or
two separate functions.
a. Set T OL = 10 12 and find the real root of the “Newton’s cubic” x3
2x
5 = 0 (the only equation
Newton ever bothered to solve with his method, in 1671 !)
b. Use it to find the three roots of the cubic (from Lab 4) x3 + x2
3x
3 = 0. Is Newton method more
efficient than Bisection?
2. Fractals from Newton Iterations. It is interesting to examine the convergence of Newton method to
find the roots for various starting points. Theory tells us that each root has a certain “region of attraction”,
such that if the initial guess lies in this region then the iteration will converge to the target root eventually.
What happens near the borders of such regions belonging to different roots? Is there a clear demarcation
between regions? Let’s explore it, for functions F (z) of a complex variable z = x+iy. Since z is representable
as a pair (x, y) of real numbers, we can visualize the regions of attraction on the plane. Consider a function
p
F (z) with several roots, e.g. F (z) = z 3 1 (its three roots are 1 and 12 ( 1 ± i 3) lying on the unit circle).
Pick a window (rectangle [a, b] ⇥ [c, d]) containing all the roots inside it, e.g. W = {(x, y) :
3
1.2 x
1.2, 1.2 y 1.2}; pick one of the roots as the “target root”, e.g. Ztgt = 1 = 1 + i ⇤ 0; pick an “escape
radius”, e.g. R = 0.2 ; pick a number of iteration to track, e.g. Nmax = 20; and finally pick a number of
points to be generated inside the window, e.g. M = 50. Subdivide [a, b] and [c, d] into M subintervals each,
thus generating a 2-dimensional grid of (M + 1) ⇥ (M + 1) points (xi , yj ), i, j = 0, ..., M . For each of these
grid points in W , we construct the following Escape-Time Algorithm for Newton iterations:
Taking a grid point z as starting value, calculate the Newton iterates
z0 = z,
zn+1 = zn
F (zn )/F 0 (zn ),
n = 0, · · · , Nmax .
If an iterate gets within R of Ztgt , then we consider it converged, and we output the coordinates xi , yj of the
starting point. If all Nmax iterates lie outside the disk of radius R centered at Ztgt , then we say the iterates
of this point escape up to “time” Nmax , and we don’t output the starting point (or print in another file for
plotting).
a. Implement this “Escape-time Algorithm” for the function F (z) = z 3
1, in a Python code. It should
accept as data: M, R, Nmax , a, b, c, d (you may hard code a, b, c, d and read the rest from a “data” file).
Set the target root in the code (so can use formulas). e.g.,
Ztgt = complex( 1.0, 0.0 )
b. To debug, use Ztgt = 1, R = 0.2, Nmax = 20, W as above, and M = 10, then M = 200. Output
the “converging” points xi , yj into a file “out1”. Plot the file from M = 200. You should see “fractal”
structures along the ±60 and 180 degree rays!
c. Play around with R, Nmax , M , describe what happens as you vary each. Record your observations.
d. Now fix the parameters: M = 400, R = 0.2, Nmax = 20, and do it for each of the roots, outputting in
separate files ’out1,txt’, ’out2.txt’ and ’out3.txt’.
e. Plot all three files together, with points or with dots. See the regions of attraction of each root by
color and observe how intricately interwinded their boundaries are! These are the Fractal sets.
f. Finally, output the plot to a PNG file (equivalent to JPG) to be included in your report.
h. Repeat a) to f) with the function F (z) = (z + 1)2 (z 2 + 0.25) and the relaxed Newton’s method
z0 = z,
(The roots are
zn+1 = zn
2F (zn )/F 0 (zn ),
n = 0, · · · , Nmax .
1, ±0.5i)
3. The Julia set associated with the complex function f (z) = z 2 +c may be depicted using the following
algorithm.
For each point, z0 , in the complex plane with
according to zn+1 = zn2 + c.
1.5 Re(z0 ) 1.5 and
1.5 Im(z0 ) 1.5, iterate
Color the pixel in an image corresponding to this region of complex plane
according to the number of iterations required for |z| to exceed some critical value, |z|max (or black if this
does not happen before a certain maximum number of iterations nmax ).
Plot the Julia set for c =
0.1 + 0.65j using |z|max = 10 and nmax = 500.
4
4
IV. Least Squares Fitting
1. A civil engineer is monitoring a leaning tower by measuring the angle it makes with a vertical line. The
engineer knows that the tower will collapse if the angle were to exceed 23o . Over the years he has recorded
the following data:
5
6
9
11
11
13.5
14
16
17
1960
1965
1970
1975
1980
1985
1990
1995
2000
Angle xi (o )
Year yi
Use the least-squares regression to fit the data to a straight line y(x) ⇡ a0 x + a1 and compute the sum of
squared error Sr , the standard error of estimate Sy/x and the coefficient of correlation r. Plot both the data
points and the linear fit on the same axes. Use the line’s equation so obtained to predict when the tower
will collapse.
Sr =
n
X1
(yk
2
a0 xk
a1 ) , Sy/x =
k=0
r
Sr
n
n
2
P
, r= q P
n x2k
P
P
xk yk ( xk )( yk )
q P
P
P
2
2
( xk ) n yk2 ( yk )
(i) Method 1. Compute the slope and intercept of the line y = a0 x + a1 using the following formulas.
a0 =
n
P
x y
Pk k
n x2k
P P
n 1
n 1
x
y
1X
1X
P k 2 k , x̄ =
xk , ȳ =
yk , a1 = ȳ
( xk )
n
n
k=0
a0 x̄
k=0
(ii) Method 2. Find the values of a0 and a1 by solving the linear system Xa = y in the least squares sense.
(X: design matrix, y: observation vector)
2
x0
6
6 x1
6
6
6 ·
6
6 ·
4
xn 1
1
3
2
y0
3
7
6
7
"
# 6 y1 7
1 7
7 a
6
7
0
7
6
7
=6 · 7
· 7
7 a1
6
7
6 · 7
· 7
5
4
5
1
yn 1
The least squares approximation to this over-determined linear system is equivalent to solving the
following system (normal equations).
(X T X)a = X T y
Define X T X = M and X T y = b and solve the 2⇥2 linear system M a = b using numpy.linalg.solve.
2. A quick way of fitting a function of the form
f (x) ⇡
a + bx
1 + cx
is to apply the least-squares method to the problem (1 + cx)f (x) ⇡ a + bx. Use this technique to fit
the world population data given here (in billions):
5
Year
Population
1000
0.340
1650
0.545
1800
0.907
1900
1.61
1950
2.56
1960
3.15
1970
3.65
1980
4.20
1990
5.30
2000
6.12
2010
6.98
2020
7.75
Plot the points and the approximated function f (x) on the same plot. Does it seem to fit well?
Determine when the world population becomes infinite!
Hint: Write a + bxi
1 and obtain a n ⇥ 3 system of equations Xa = y and
cxi yi = yi , i = 0, 1, 2, · · · , n
solve it in the least square sense [(X X)a = X y]
2
3
2
3
3
1
x0
x 0 y0 2
y0
6
7 a
6
7
6 1
7
x1
x 1 y1 7
y1
7 6
6
76
6
7
=6 .
6 .
76
7
b 7
.
.
4
5
6 ..
6 ..
7
..
.. 7
4
5 c
4
5
1 xn 1
x n 1 yn 1
yn 1
T
5
T
V. Numerical Methods for Integration and ODEs
1. Implement the following numerical methods to approximate the definite integral
Z 1
2
e x dx.
0
a. Midpoint rule:
Z b
a
f (x)dx ⇡ h
b. Trapezoid rule:
Z b
a
N
X
f (xi ), h =
i=1
b
a
N
, xi = a +
"
#
N
X1
1
1
f (x) dx ⇡ h f (x0 ) +
f (xi ) + f (xN ) ,
2
2
i=1
h=
2i
1
2
b
a
N
h.
,
xi = a + ih
c. Simpson’s rule:
2
3
N
N
X1
X2
h4
b a
f (x)dx ⇡
f (x0 ) + 4
f (xi ) + 2
f (xi ) + f (xN )5 , h =
, xi = a + ih.
3
N
a
i=1,3,5
i=2,4,6
Z b
d. 2-point Gauss quadrature rule:
Z b
a
N
f (x)dx ⇡
hX
1
[f (xl ) + f (xr )] , xl =
2 i=1
2
h=
b
✓
a
N
6
xi + xi 1
, xi = a + ih
h
p
3
◆
, xr =
1
2
✓
h
xi + xi 1 + p
3
◆
p
Comparing against the true solution 12 ⇡erf(1), which method gives the best approximation to ⇡?
How accurate did you get? Display your results (error of approximation for different values of N ) both
graphically and in tabular form.
2. Length and width measurements (in feet) of an irregularly shaped land is shown in the table below.
Suppose you want to lay down a uniform 3-inch covering of pine bark mulch over the entire plot. Home
improvement store sells bags containing 3 cubic feet of bark mulch, how many bags do you need? Use
the Simpson’s rule to approximate area.
x
0
1
2
3
4
5
y
8
10 11 15 16
6
7
8
9
10
11 12
13
14
15
16
16 16 16
15.5
15.5
15.5
15 15
15 14.5
14.5
14
x
17 18
19
20 21 22
23 24 25
26
y
14 14
13.5
13 13 13
12 11 19
6
3. The logistic model is used to simulate population as in
✓
◆
dp
p
= krm 1
p
dt
pmax
where p(t) is the population at time t, krm is the maximum growth rate under unlimited conditions, and pmax is the carrying capacity. Simulate the world’s population from 1950 to 2020 using
the Euler method. Employ the following initial conditions and parameter values for your simulation:
p0 (in 1950) = 2555 million people, krm = 0.026/yr, and pmax = 12, 000 million people. Have the function generate output corresponding to the dates for the following measured population data. Develop
a plot of your simulation along with the data.
t
1950
1960
1970
1980
1990
2000
2010
2020
p
2555
3040
3708
4454
5276
6079
6922
7753
7
c
MA305.01 – Khanal
Fall 2019
GUIDE on Preparing Final Report
The Final Report should be as self-contained as possible, readable, to the point, understandable by an educated non-expert. It should be typed, preferably using LaTeX
(a text-formatting tool) and organized as follows.
0. Title page
1. Introduction
2. Problem description
3. Methods
4. Solution
5. Discussion/Conclusions
6. References/Appendix
7. Acknowledgement
Here are some general guidelines for each part. These are guidelines, try to follow them
to the extent that’s reasonable for each project. Some parts may apply more or less to
each particular project, ... use your judgement!
0. Title page
MA305.01 – Final Project
Project Title
Your Name(s), Date
1. Introduction
Write a very brief description (background/significance) of what the project is about.
2. Problem Statement and Assumptions
State any assumptions made for the formulation of the model. State fully and
precisely the mathematical problem, Explain meaning of all symbols used, Make
clear what is given and what we are looking for.
3. Method
Begin with naming or characterizing the method/approach to be used, perhaps explain the basic idea behind it, to what type of problems it applies, under what
conditions, what it achieves, what are its main features, advantages, disadvantages.
Justify why it is applicable to this problem, stating clearly any assumptions you
need to make about the problem for the method to apply. Name some other methods/approaches one could use, and if/why your method may be preferable.
4. Solutions/Results (Presentation of your solution and results)
Describe your implementation of the method(s) for this specific problem, any special
features, parallelization strategy, choices of any parameters, stopping criteria, etc.
Present the results in words and plots (annotate by hand if necessary), explain what
they mean. Include your code in an Appendix.
5. Discussion/Conclusions
Interpret your solution physically, what we learn from it, comment on strengths/weaknesses
of the solution method, any nice features you want to brag about, possible ways to
improve it (e.g. how to make it more accurate, more efficient), as appropriate.
6. References/Appendix
List the materials used in the project. e.g., books, papers, web resources, codes, etc.
7. Acknowledgement
If there are any acknowledgments to research support etc., they go here.
MA305, Fall 2019
Embry-Riddle Aeronautical
University
E
Title of the Project
Your Name(s)
December 3, 2019
Abstract
Contents
1 Introduction
3
2 Problem Statement
2.1 Part I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Part II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Part III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3
3
3
3 Method/Analysis
3
4 Solutions/Results
4.1 A subsection . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 A subsubsection . . . . . . . . . . . . . . . . . . . . .
4.1.2 A further subdivision . . . . . . . . . . . . . . . . . .
3
4
4
4
5 Discussion/Conclusions
4
A Python Codes
5
2
1
Introduction
Write a brief description (background/significance) of what the project is
about.
2
Problem Statement
State fully and precisely the mathematical problem. Explain meaning of all
symbols used. Make clear what is given and what we are looking for.
2.1
Part I
Text introducing this subsection.
2.2
Part II
Text introducing this subsection.
2.3
Part III
Text introducing this subsection.
3
Method/Analysis
Begin with naming or characterizing the method/approach to be used, perhaps explain the basic idea behind it, to what type of problems it applies,
under what conditions, what it achieves, what are its main features, advantages, disadvantages. Justify why it is applicable to this problem, stating
clearly any assumptions you need to make about the problem for the method
to apply. Name some other methods/approaches one could use, and if/why
your method may be preferable.
4
Solutions/Results
This section contains the presentation of your solution and results. Describe your implementation of the method(s) for this specific problem, any
special features, numerical methods implementation strategy, choices of any
parameters, stopping criteria, etc. Present the results in words and plots
(annotate by hand if necessary), explain what they mean. Include your code
in an Appendix.
3
4.1
A subsection
Text introducing this subsection.
4.1.1
A subsubsection
Text introducing this subsubsection.
4.1.2
A further subdivision
Text introducing this subsubsection.
5
Discussion/Conclusions
Interpret your solution physically, what we learn from it, comment on strengths
and weaknesses of the solution method, any nice features you want to brag
about, possible ways to improve it (e.g. how to make it more accurate, more
efficient), as appropriate.
References
[1] Heath, Michael T., Scientific Computing: An Introductory Survey,
McGraw Hill, 2002.
4
A
Python Codes
Text introducing this appendix. Subsections and further divisions can also
be used in appendices.
1 #! / u s r / b i n / env
2
python3
import math
3
4
”””
5 =================================================================
MA305
cw #: your name
date
Purpose : Find t h e number o f t r a i l i n g z e r o s i n f a c t o r i a l ( n ) .
8 =================================================================
6
7
9
”””
10
11
12
n=i n p u t ( ’ Enter a p o s i t i v e i n t e g e r : ’ )
n=i n t ( n )
13
count=0
f o r i i n r a n g e ( 1 , n+1) :
16
count += n //5⇤⇤ i #pow ( 5 , i )
14
15
17
18
19
p r i n t ( ’ Number o f t r a i l i n g z e r o s i n f a c t o r i a l ( ’ , n , ’ ) : ’ , count )
p r i n t ( ’ F a c t o r i a l ( ’ , n , ’ )= ’ , math . f a c t o r i a l ( n ) )
5