1. shopping
2. resources
3. coupons

A Solution Manual and Notes for:
Numerical Methods in Engineering with Python
by Jaan Kiusalaas
John L. Weatherwax∗
October 6, 2014
∗
[email protected]
1
Chapter 2 (Systems of Linear Algebraic Equations)
Notes on Section 2.1
Notes on the Choleski decomposition
In the derivation of the Choleski decomposition we wanted to compute (LLT )ij which we
can do with the following algebraic steps
(LLT )ij =
n
X
Lik (LT )kj =
n
X
Lik Ljk .
k=1
k=1
Because L is lower triangular we know that Ljk = 0 when k > j and the above becomes
j
X
Lik Ljk ,
k=1
which is the expression given in the book for (LLT )ij .
Problem Set 2.1
Problem 1
To evaluate the determinant of these matrices I will perform an LU decomposition and then
take the product of the diagonal elements of the matrix U. To do this I will use the python
code LUdecomp. This problem is worked in the python code c 2 1 p 1.py. When we run
that code we get
Part
Part
Part
Part
(a):
(b):
(c):
(d):
determinant=
determinant=
determinant=
determinant=
-0.000000,
0.647537,
4.000000,
-0.000000,
determinant/mean(Aij)=
determinant/mean(Aij)=
determinant/mean(Aij)=
determinant/mean(Aij)=
-0.000000
0.267209
3.600000
-0.000000
This indicates that the matrices in Part (a) and (d) are singular (or very poorly ill-conditioned)
and that the matrices from Part (b) and (c) are well-conditioned.
Problem 2
To compute the matrix A given L and U we would just compute the matrix product LU. To
determine the determinant of A we need to take the product of the elements on the diagonal
of the U matrix.
2
Problem 3
Given the LU decomposition we first need to solve Ly = b by forward substitution and then
Ux = y by back substitution. The system Ly = b is


 

1
0
0
y1
1
 3/2
1
0   y2  =  −1  .
1/2 11/13 1
y3
2
First we have y1 = 1. Next for y2 we have
3
5
3
y2 = −1 − y1 = −1 − = − .
2
2
2
Finally for y3 we have
5
47
11
1 11
1
−
=
.
y3 = 2 − y1 − y2 = 2 − −
2
13
2 13
2
13
Next the system Ux = y is


 

2 −3
−1
x1
1
 0 13/2 −7/2   x2  =  −5/2  .
0
0
32/13
x3
47/13
47
.
32
Next for x2 we find
2
5 32
5 32 47
159
2
x2 =
− − x3 =
− −
=−
.
13
2 13
13
2 13 32
169
Finally we compute
1746
.
x1 = −
1061
Thus x3 =
Problem 5
We will use the routine gaussElimin to solve for x for one of the right-hand-sides for this
problem. This problem is worked in the python code c 2 1 p 5.py. When we run that code
we get
Solution=
[ 0.4375
0.25
-0.125
-0.0625]
Problem 6
We will use the routine gaussElimin to solve for x for this problem. Because of the way
this code is written we have to perform “pivoting” by hand (which we do in the code for
this problem). This problem is worked in the python code c 2 1 p 6.py. When we run that
code we get
Solution=
[ 2. -2.
1.
1. -1.]
3
Problem 7
Part (a): For this part we will use the routine LUdecomp to find L and U in this problem.
This problem is worked in the python code c 2 1 p 7.py. When we run that code we get
A after LUdecomp (Doolittle)=
[[ 4.
-1.
0.
]
[-0.25
3.75
-1.
]
[ 0.
-0.26666667 3.73333333]]
This means that our matrix L and U such that A = LU are given by


1
0
0
1
0 
L =  −0.25
0
−0.26666667 1


4 −1
0
.
−1
U =  0 3.75
0 0 3.73333333
Part (b): This part is worked in the routine choleski. Running the above code gives
A after Choleski=
[[ 2.
0.
[-0.5
1.93649167
[ 0.
-0.51639778
0.
]
0.
]
1.93218357]]
Note that this is the matrix L in the Choleski decomposition A = LLT .
Problem 8
This problem is worked in the python code c 2 1 p 8.py and uses the routines LUdecomp
and LUsolve to perform Doolittle’s LU decomposition on A. When we run that code we get
the following
Solution =
[ 1.
2.
3.]
Problem 9
This problem is worked in the python code c 2 1 p 9.py and uses the routines LUdecomp
and LUsolve to perform Doolittle’s LU decomposition on A. When we run that code we get
the following
4
Solution =
[ 1.
1.
1.]
Problem 10
This problem is worked in the python code c 2 1 p 10.py and uses the routines LUdecomp
and LUsolve to perform Doolittle’s LU decomposition on A. When we run that code we get
the following
Solution to Ax=B[:,0] is
Solution to Ax=B[:,1] is
[-3.75
[ 1.75
1.66666667 3.5
-0.66666667 -1.5
]
]
Problem 11
This problem is worked in the python code c 2 1 p 11.py and uses the routines choleski.
When we run that code we get the following
Solution=
[ 0.5 -1.
1.5]
Problem 12
This problem is worked in the python code c 2 1 p 12.py and uses the routines LUdecomp
and LUsolve to perform Doolittle’s LU decomposition on A. When we run that code we get
the following
Solution=
[ 0.98
-0.22125
1.0875 ]
Problem 13
For a diagonal matrix to be positive definite means that αi > 0 for all i. Then for A = LLT
√
and to have A diagonal means that Lij = 0 for all i 6= j and Lii = αi for the diagonal
elements.
Problem 14
This problem is worked by calling the routine gaussEliminMultipleRHS from the python
code c 2 1 p 14.py When we run that code we get the following
5
Solution=
[ 1. 2.
[ 1. 2.
[[ 1.
2.]
3.]]
1.
1.]
Problem 15
This problem is worked in the python code c 2 1 p 15.py. When we run that code we get
the following
n=
n=
n=
n=
n=
n=
n=
3;
4;
5;
6;
7;
8;
9;
error=
error=
error=
error=
error=
error=
error=
0.000000
0.000000
0.000000
0.000000
0.000000
0.000001
0.000036
When n = 9 the error is greater than 10−6 and the program stops.
Problem 16
This is already implemented in the choleski code in a function called choleskiSol.
Problem 17
This problem is equivalent to a linear system y = Ax where
y=
y1 y2 y3 y4
the four given yi values, A is given by

1
 1
A=
 1
1
x1
x2
x3
x4
x21
x22
x23
x24
T
,

x31
x32 
,
x33 
x34
T
and the unknown x is the vector of coefficients x = a0 a1 a2 a3 . This problem is
worked in the python code c 2 1 p 17.py. When we run that code we get the following
Solution=
[10.0, 34.0, -9.0, 0.0]
6
which means that the best fit polynomial is given by
y = 10 + 34x − 9x2 .
Problem 18
See the previous problem for the formulation used. This problem is worked in the python
code c 2 1 p 18.py. When we run that code we get the following
Solution=
[-1.0, 2.6833, -0.875, 0.2167, -0.025]
which means that the best fit polynomial is given by
y = −1 + 2.6833x − 0.875x2 + 0.2167x3 − 0.025x4 .
Problem 19
Our function y(x) takes the form
y(x) = a0 + a1 x + a2 x2 + a3 x3 + a4 x4 .
The curvature of a curve y = y(x) is proportional to the second derivative of y(x) thus to
make the curvature vanish we can make y ′′(x) vanish. Computing the first two derivatives
of y(x) we have
y ′(x) = a1 + 2a2 x + 3a3 x2 + 4a4 x3
y ′′(x) = 2a2 + 6a3 x + 12a4 x2 .
This means that we want the coefficients of our polynomial y(x) to satisfy the conditions
y(0) = 1
y(0.75) = −0.25
y(1) = 1
y ′′ (0) = 0
y ′′ (1) = 0 .
From these the linear system Ax = b we need to solve has the vector b given by
b=
a vector of unknowns x given by
x=
1 −0.25 1 0 0
T
,
a0 a1 a2 a3 a4
T
,
7
and a matrix A given by



A=


1 0
02
03
04
1 0.75 0.752 0.753 0.754
1 1
12
13
14
0 0
2
6(0) 12(02 )
0 0
2
6(1) 12(12 )



.


Using these this problem is worked in the python code c 2 1 p 19.py. When we run that
code we get the following
Solution=
[1.0, -5.614, 0.0, 11.2281, -5.614]
which means that the best fit polynomial that satisfies all of the required conditions is given
by
y = 1.0 − 5.614x + 11.2281x3 − 5.614x4 .
Problem 20
This problem is worked in the python code c 2 1 p 20.py. When we run that code we get
the following
|A|= -0.225797, |A|/mean(abs(A))=
1.495072
b-Ax= [ -1.77635684e-15
8.88178420e-16
0.00000000e+00
0.00000000e+00]
Using the metric that an ill-conditioned matrix is one with a determinant that is much smaller
than the elements of the matrix I don’t see this matrix as that ill-conditioned. Computing
Ax and comparing it with b shows that the two vectors are really quite close.
8
Chapter 3 (Interpolation and Curve Fitting)
Problem Set 3.1
Problem 1
Part (a): We could use the routine neville to do this part but we will instead do the
calculations by hand. Do do this we use the formula in the book for recursively computing Pk [xi , xi+1 , . . . , xi+k ] in terms of Pk−1 [xi , xi+1 , . . . , xi+k−1 ] and Pk−1 [xi+1 , xi+2 , . . . , xi+k−1 ].
We start with
P0 [x0 ] = y0 = −5.76
P0 [x1 ] = y1 = −5.61
P0 [x2 ] = y2 = −3.69 .
Next we compute
(x − x1 )P0 [x0 ] + (x0 − x)P0 [x1 ]
x0 − x1
(x − 0.3)(−5.76) + (−1.2 − x)(−5.61)
=
.
−1.2 − 0.3
P1 [x0 , x1 ] =
If x = 0 this becomes P1 [x0 , x1 ] = −5.64. Next we compute
(x − x2 )P0 [x1 ] + (x1 − x)P0 [x2 ]
x1 − x2
(x − 1.1)(−5.61) + (0.3 − x)(−3.69)
.
=
0.3 − 1.1
P1 [x1 , x2 ] =
If x = 0 this becomes P1 [x1 , x2 ] = −6.33. Next we compute
(x − x2 )P1 [x0 , x1 ] + (x0 − x)P1 [x1 , x2 ]
x0 − x2
(x − 1.1)(−5.64) + (−1.2 − x)(−6.33)
.
=
−1.2 − 1.1
P2 [x0 , x1 , x2 ] =
If x = 0 this becomes P2 [x0 , x1 , x2 ] = −6. We can check that our result is correct by using
the neville code. We do this in the code c 3 1 p 1.py where we also get the value -6.0.
Part (b): In the three point case our Lagrange interpolating polynomial is given by
P2 (x) = y0 l0 (x) + y1 l1 (x) + y2 l2 (x)
(x − x1 )(x − x2 )
(x − x0 )(x − x2 )
(x − x0 )(x − x1 )
= y0
+ y1
+ y2
(x0 − x1 )(x0 − x2 )
(x1 − x0 )(x1 − x2 )
(x2 − x0 )(x2 − x1 )
(x + 1.2)(x − 1.1)
(x + 1.2)(x − 0.3)
(x − 0.3)(x − 1.1)
− 5.61
− 3.69
= −5.76
(−1.2 − 0.3)(−1.2 − 1.1)
(0.3 + 1.2)(0.3 − 1.1)
(1.1 + 1.2)(1.1 − 0.3)
= −1.669565(x − 0.3)(x − 1.1) + 4.675(x + 1.2)(x − 1.1) − 2.005435(x + 1.2)(x − 0.3) .
If we let x = 0 we get P2 (0) = −6.
9
Problem 2
For this problem we will use inverse interpolation to model the function x = x(y). Using
this expression and taking y = 0 gives us an approximate root.
Part (a): The three closes values to y = 0 are the pairs (in (y, x) notation)
(−1.5727, 3) , (−0.4112, 2.5) , (0.8509, 2) .
With these points the Lagrangian interpolating polynomial is given by
(y − y0 )(y − y2 )
(y − y0 )(y − y1 )
(y − y1 )(y − y2 )
+ x1
x3
(y0 − y1 )(y0 − y2 )
(y1 − y0 )(y1 − y2 ) (y2 − y0 )(y2 − y1 )
(y + 0.4112)(y − 0.8509)
(y + 1.5727)(y − 0.8509)
=3
+ 2.5
(−1.5727 + 0.4112)(−1.5727 − 0.8509)
(−0.4112 + 1.5727)(−0.4112 − 0.8509)
(y + 1.5727)(y + 0.4112)
+2
(0.8509 + 1.5727)(0.8509 + 0.4112)
= 1.0657150(y + 0.4112)(y − 0.8509) − 1.7054030(y + 1.5727)(y − 0.8509)
+ 0.6538457(y + 1.5727)(y + 0.4112) .
x(y) = x0
If we take y = 0 in the above we find x(y = 0) = 2.332143.
Part (a): To use four points we would need to include the point (1.9047, 2) in the above
model. This would require us to add another factor in the numerator and the denominator
of the Lagrangian basis functions we already have and one more basis function all together.
Problem 3
See the code c 3 1 p 3.py where we get the value 2.55710037008.
Problem 4
See the code c 3 1 p 4.py where we get the value 3.09553870532.
Problem 5
See the code c 3 1 p 5.py where we computed the value of y at each of the given x values
using three methods: a newton interpolating polynomial, Neville’s method, and cubic splines.
When we run that code we get
x=
x=
0.785398; newton=
1.570796; newton=
1.336495, neville=
2.214944, neville=
10
1.336495, spline=
2.214944, spline=
1.333127
2.205415
Problem 6
See the code c 3 1 p 6.py. For this problem we modified the in-place algorithm in the book
for computing the Newton divided difference table to build the full table. The code then
prints this table. When we run the above code we get
[[ -1.
[ 2.
[ 59.
[ 4.
[ 24.
[-53.
0.
1.
10.
5.
5.
26.
0.
0.
3.
-2.
2.
-5.
0.
0.
0.
1.
1.
1.
0.
0.
0.
0.
0.
0.
0.]
0.]
0.]
0.]
0.]
0.]]
Based on the above and the correspondence between the diagonal elements in the above grid
and the coefficients an in Newton’s polynomial representation, we see that a0 , a1 , a2 , and a3
are non-zero but an for n ≥ 4 are. Thus the given data is fit by a cubic polynomial.
Problem 7
See the code c 3 1 p 7.py. When we run that we get
Coefficients from newtonPoly.coeffts=
[0 1 1 0 0]
Thus the polynomial that fits this data is given by
P2 (x) = (x + 3) + (x + 3)(x − 2) .
Problem 8
To use Neville’s method to evaluate the equation of the quadratic that passes through the
points we want to evaluate the expression P2 [x0 , x1 , x2 ] as a function of x i.e. the polynomial
itself and not the value of the polynomial evaluated at a specific numerical x value. Note
that the routine neville computes the value of the interpolating polynomial at a specific x
value and not the coefficients of the polynomial. To get the full polynomial we start with
P0 [x0 ] = y0 = 17
P0 [x1 ] = y1 = −7
P0 [x2 ] = y2 = −15 .
11
Next for P1 [x0 , x1 ] we compute
(x − x1 )P0 [x0 ] + (x0 − x)P0 [x1 ]
x0 − x1
17(x − 1) + 7(x + 1)
17(x − 1) + (−1 − x)(−7)
=
= −12x + 5 .
=
−1 − 1
−2
P1 [x0 , x1 ] =
Next for P1 [x1 , x2 ] we compute
(x − x2 )P0 [x1 ] + (x1 − x)P0 [x2 ]
x1 − x2
(x − 3)(−7) + (1 − x)(−15)
= −4x − 3 .
=
1−3
P1 [x1 , x2 ] =
Finally we compute
(x − x2 )P1 [x0 , x1 ] + (x0 − x)P1 [x1 , x2 ]
x0 − x2
(x − 3)(−12x + 5) + (−1 − x)(−4x − 3)
= 2x2 − 12x + 3 .
=
−1 − 3
P2 [x0 , x1 , x2 ] =
Lets check that this quadratic passes through the given points. Evaluating this expression
at x = −1, x = +1, and x = +3 we find
P2 [x0 , x1 , x2 ](−1) = 17
P2 [x0 , x1 , x2 ](+1) = −7
P2 [x0 , x1 , x2 ](+3) = −15 ,
as they should
Problem 9
The Lagrange interpolating polynomial is
(h − h1 )(h − h2 )
(h − h0 )(h − h2 )
(h − h0 )(h − h1 )
+ ρ1
+ ρ2
(h0 − h1 )(h0 − h2 )
(h1 − h0 )(h1 − h2 )
(h2 − h0 )(h2 − h1 )
(h − 3)(h − 6)
(h − 0)(h − 6)
(h − 0)(h − 3)
= 1.225
+ 0.905
+ 0.652
(0 − 3)(0 − 6)
(3 − 0)(3 − 6)
(6 − 0)(6 − 3)
= 0.06805556(h − 3)(h − 6) − 0.10055556h(h − 6) + 0.03622222h(h − 3) ,
ρ(h) = ρ0
which could be simplified further if needed.
Problem 10
When we use only three points our cubic spline is made up of two pieces f0,1 (x) and f1,2 (x)
where we can use Eq. 3.10 to give fi,i+1 (x) in terms of xi , yi , and ki ’s. Using the curvatures
function in the cubicSpline module we compute values for ki given xi and yi (see the code
c 3 1 p 10.py) to get
12
k=
[ 0.
-4.5
0. ]
Using these we find our two splines given by
y0 (x − x1 ) − y1 (x − x0 )
k1 (x − x0 )3
− (x − x0 )(x0 − x1 ) +
f0,1 (x) = −
6 x0 − x1
x0 − x1
3
3
11
= (−x3 + x) + 2x = − x3 + x ,
4
4
4
and
k1 (x − x2 )3
y1 (x − x2 ) − y2 (x − x1 )
f1,2 (x) =
− (x − x2 )(x1 − x2 ) +
6 x1 − x2
x1 − x2
3
2(x − 2) − (x − 1)
3 (x − 2)
+ (x − 2) +
=−
4
(−1)
(−1)
3
11
= (x − 2)3 − (x − 2) + (x − 1) .
4
4
Given these expressions we compute
11
1
9
′
′
so f0,1
(1) =
f0,1
(x) = − x2 +
4
4
2
9
9
′′
′′
f0,1 (x) = − x so f0,1 (1) = − ,
2
2
and
11
1
9
′
′
+ 1 so f0,1
(1) =
f1,2
(x) = (x − 2)2 −
4
4
2
9
9
′′
′′
f1,2 (x) = (x − 2) so f0,1 (1) = − ,
2
2
showing that the cubics have the same first and second derivatives.
Problem 11
See the code c 3 1 p 11.py. When we run that we get
x=
3.400000; spline=
10.254857
Problem 12
See the code c 3 1 p 12.py. When we run that we get
approximate zero (inverse interpolation)=
13
0.721157
Problem 13
Example 3.6 has four data points (so n = 3) of unit spacing h = 1 between the xi . From
Eq 3.12 we have n − 1 equations for the n + 1 unknowns ki for i = 0, 1, 2, . . . , n − 1, n. We
thus need two additional conditions to make this a solvable system which for this problem
will be k0 = k1 and k2 = kn−1 = kn = k3 . Let i = 1 and i = n − 1 = 2 in Eq. 3.12 to get
k0 + 4k1 + k2 = 6(y0 − 2y1 + y2 )
k1 + 4k2 + k3 = 6(y1 − 2y2 + y3 ) .
Since we know that k0 = k1 and k2 = k3 we can write these two equations as
5k1 + k2 = 6(y0 − 2y1 + y2 )
k1 + 5k2 = 6(y1 − 2y2 + y3 ) .
This is thus a system of two equations with two unknowns which we can solve for to find
(k1 , k2 ). When we do that (for the given values of yi ) we get k1 = −0.625 and k2 = 0.125.
The vector k of curvatures needed in the python implementation is given by
k = [ -0.625, -0.625, 0.125, 0.125 ]
We can pass a vector like this to the evalSpline routine to evaluate our spline at auxiliary
points if needed.
Problem 14
See the code c 3 1 p 14.py. When we run that we get
x=
x=
x=
1.100000; neville=
1.200000; neville=
1.300000; neville=
1.326194
1.393758
1.469308
Note we could have used the python command raw input to have the user enter points to
compute the interpolant at.
Problem 15
See the code c 3 1 p 15.py. When we run that we get
T= 200.000000; neville=
T= 400.000000; neville=
0.993335, spline=
0.985982, spline=
14
0.991707
1.088293
Problem 16
See the code c 3 1 p 16.py. When we run that we get
T=
0.460000; neville=
4.878536, spline=
3.843401
Problem 17
See the code c 3 1 p 17.py. When we run that we get
Re=
Re=
Re=
Re=
5.000;
50.000;
500.000;
5000.000;
spline=
spline=
spline=
spline=
6.901592
1.590835
0.557434
0.386822
Problem 18
See the code c 3 1 p 18.py. When we run that we get
Re=
Re=
Re=
Re=
5.000;
50.000;
500.000;
5000.000;
NN
NN
NN
NN
Neville=
Neville=
Neville=
Neville=
6.932565
1.596500
0.552888
0.372227
Problem 19
See the code c 3 1 p 19.py. When we run that we get
T=
T=
T=
T=
10.000000;
30.000000;
60.000000;
90.000000;
neville=
neville=
neville=
neville=
1.620692,
0.842478,
0.457169,
0.333423,
spline=
spline=
spline=
spline=
Problem 20
See the code c 3 1 p 20.py. When we run that we get
15
1.474981
0.870148
0.454081
0.319453
h=
10.500000; neville=
0.317776, spline=
0.309018
Problem Set 3.2
Problem 1
The equation for a (this is the books Eq. 3.19) is y¯ − x¯b. This means that
y¯ = a + b¯
x,
which means that the point (¯
x, y¯) is on the line.
Problem 2 (simple least-squares)
For this problem we can follow Example 3.10 from the text and using polyFit. We implement
this in the python code c 3 2 p 2.py which gives the following
coeff= [-0.02 1. ]
std= 0.0316227766017
This means that the least-square (best fit line) is
y = −0.02 + x ,
and that the spread of the data about this line is given by σ = 0.0316227766017.
Problem 3 (ordinary least-squares)
In this problem we have several measurements (three) of the strain when a given stress is
applied. As no machine is stated to have better measurement accuracy than the others We
can attempt to solve this problem using ordinary least-squares (not weighted least-squares).
To do that we repeat the four stress measurements three times (once for each of the three
strain tests) and then determine the coefficient a in the model
strain = a stress + b .
We would expect b ≈ 0 (since when stress = 0 we expect to have strain = 0) and should
really model this data in that way but since the book provides no background on how to
model data without an intercept we will simple model the data as a first order polynomial
using polyFit. We implement all of this in the python code c 3 2 p 3.py which gives the
following output
16
coeff= [ -1.66666667e-06
std= 0.000110819673344
1.48502415e-11]
The second number in the coeff output or 1.48502415 10−11 is the estimate of the modulus
of elasticity.
Problem 4 (weighted least-squares)
Note that many of the comments in Problem 3 hold for this problem as well. In this case
since the third machine is less accurate than the other two we should model it as such. That
is we want to find coefficients in our model a0 , a1 , . . . , am that minimize
S(a0 , a1 , . . . , am ) =
n
X
i=0
=
wi2 (yi − f (xi ))2
X
(yi − f (xi ))2 +
machine 1
(1)
X
(yi − f (xi ))2 +
machine 2
X
machine 3
2
1
(yi − f (xi ))2 .
2
(2)
Note that in problem 3 above we took wi = 1 in Equation 1 i.e. did not have any weights to
worry about. In this problem this expression take the form of Equation 2 where the third
machine has one-half of the weight of the other two. We can solve such a problem using
weighted least-squares. We implement all of this in the python code c 3 2 p 4.py which
gives the following output
coeff=
[
9.18918919e-06
1.44402664e-11]
The second number in the coeff output or 1.44402664 10−11 is the estimate of the modulus
of elasticity.
Problem 5
This problem is worked in the python code c 3 2 p 5.py which gives the following output
coeff= [ 3.00559091 -0.43838182]
std= 0.0774987715771
This means the linear fit found is
y = 3.00559091 − 0.43838182x ,
17
Problem 6
This problem is worked in the python code c 3 2 p 6.py which gives the following output
coeff= [ 1.86191098e+01
std= 0.559659297904
-5.89632543e-03]
Which means that our least-square line for this data is given by
φ = 18.6191098 − 0.00589632543M .
Problem 7
This problem is worked in the python code c 3 2 p 7.py which gives the following output
coeff= [ 0.99889524 -0.09344731 0.00276321]
p(x),p’(x),p’’(x) at x=10.5= (0.32234242971039695, (-0.035419894805543056+0j), (0.005526420189145401+0j))
Problem 8
This problem is worked in the python code c 3 2 p 8.py which gives the following output
coeff= [
T=
10;
T=
30;
T=
60;
T=
90;
1.79570895e+00 -3.93212797e-02
mu_k(T)=
1.434507
mu_k(T)=
0.888944
mu_k(T)=
0.436570
mu_k(T)=
0.301553
3.28569164e-04
-8.45886166e-07]
We can compare these predictions to the ones made in Problem 19 on Page 15 where we see
that they are quite similar.
Problem 9
This problem is worked in the python code c 3 2 p 9.py which gives the following output
degree=
degree=
1; sigma=
2; sigma=
2.243564 with coeff= [-6.18989525 9.43854354]
0.812928 with coeff= [ 4.40567377 -1.06889613 2.10811822]
If we select the model (the polynomial degree) based on the smaller value of σ we would
select the second degree polynomial. To really select the polynomial we should consider a
model selection technique as discussed in [1].
18
Problem 10
This problem is worked in the python code c 3 2 p 10.py which gives the following output
degree=
degree=
degree=
degree=
degree=
degree=
degree=
1;
2;
3;
4;
5;
6;
7;
sigma=
sigma=
sigma=
sigma=
sigma=
sigma=
sigma=
2.855202
2.768407
2.265845
2.234047
2.495617
17.594028
3.747578
Results like these might lead us to select the fourth order polynomial (it has the smallest
value of σ). As mentioned in the text, however, selecting a model based on its in-sample
performance can be misleading (see [1]).
Problem 11
This problem is worked in the python code c 3 2 p 11.py which gives the following output
degree=
2; sigma=
0.000418; coeff= [
1.05258944e+00 -6.80740001e-04
2.30985609e-07]
Which means that the best fit quadratic is given by
k = 1.05258944 − 0.000680740001T + 0.000000230985609T 2 .
Problem 12
The solution to a similar problem is worked in this chapter and this section of the text
entitled “Fitting exponential functions”. Following the logic there we recall that ri and Ri
are defined by
ri ≡ yi − axbi
Ri ≡ ln(yi ) − (ln(a) + b ln(xi )) .
Then since ln(axbi ) = ln(ri − yi ) we have
Ri = ln(yi ) − ln(yi − ri ) = ln
ri
ri
≈ ,
= − ln 1 −
yi
yi
yi
y i − ri
= − ln
y i − ri
yi
as we were to show. In the above we have used the approximation that ln(1+x) = x+O(x2 ).
19
Problem 13
The error sum of squares we want to minimize is given by
πx 2
πx X
i
i
− b cos
.
S(a, b) =
yi − a sin
2
2
i
To minimize S(a, b) as a function of a and b we need to take the derivatives of S with respect
to a and b and set the results equal to zero. Taking these derivatives we get
πx πx πx X
∂S
i
i
i
− b cos
− sin
=0
=2
yi − a sin
∂a
2
2
2
i
πx πx πx X
∂S
i
i
i
− b cos
− cos
= 0.
=2
yi − a sin
∂b
2
2
2
i
This as a system of 2 × 2 equations for the unknowns a and b that we can write as
πx X
πx πx πx 2
X
X
i
i
i
i
cos
=
yi sin
+b
sin
a
sin
2
2
2
2
i
i
i
πx πx πx 2 X
πx X
X
i
i
i
i
a
sin
cos
+b
cos
.
=
yi cos
2
2
2
2
i
i
i
Given data (xi , yi ) for i = 0, 1, . . . , n we can form the above system and then solve it for a
and b. This problem is worked in the python code c 3 2 p 13.py which gives the following
output
a=
3.038491; b=
-2.049560
Problem 14 (a transformation to simple linear regression)
To use least-squares regression to fit this model we take the logarithm of both sides to get
ln(y) = ln(a) + b ln(x). We then fit a simple linear regression model to the data points
(ln(x), ln(y)). Once we have the parameters from this linear fit (say aˆ and ˆb) the parameters
in the desired model are then a = eaˆ and b = ˆb. This problem is worked in the python code
c 3 2 p 14.py which gives the following output
coeff_hat= [ 0.53247348 1.88175018]
coeff= 1.70313978716 1.88175018137
which means that our estimated model is
y = 1.70313978716x1.88175018137 .
20
Problem 15 (another transformation to simple linear regression)
y
x
= aebx , and then fit this model by taking logarithms of
y = ln(a) + bx .
ln
x
If we call the estimates of ln(a) and b when we fit the points x, ln xy the aˆ and ˆb then the
parameters we want are given by a = eaˆ and b = ˆb. This problem is worked in the python
code c 3 2 p 15.py which gives the following output
First divide by x to get the model
both sides
coeff_hat= [ 1.07196796 -1.98387741]
coeff= 2.92112249311 -1.98387740743
sigma= 0.00599685706921
which means that our estimated model is
y = 2.92112249311xe−1.98387740743x ,
and that σ has the given value.
Problem 16
To use least-squares regression to fit this model we take the logarithm of both sides to get
ln(γ) = ln(a) − bt. We then fit a simple linear regression model to the data points (t, ln(γ)).
Once we have the parameters from this linear fit (say a
ˆ and ˆb) the parameters in the desired
a
ˆ
model are then a = e and b = −ˆb. Once we have these parameters the estimate of the
half-life T1/2 is given by
ln(2)
T1/2 =
.
b
This problem is worked in the python code c 3 2 p 16.py which gives the following output
coeff_hat= [-0.00158547 -0.00863955]
coeff= 0.998415781283 0.00863954970145
T_{1/2}= 80.229549515
This means that our model is estimated by
γ(t) = 0.998415781283e−0.00863954970145t ,
with the given half-life (in years).
21
80
60
f(x)
40
20
0
−20
−40
−3
−2
−1
0
1
2
3
4
5
6
x
Figure 1: A plot of the function f (x) for Problem 2.
Chapter 4 (Roots of Equations)
Problem Set 4.1
Problem 1
For this problem we will look for a root of the function f (x) = x3 − 75 using the NewtonRaphson method. This method implements the iterative scheme
3
xi − 75
f (xi )
,
= xi −
xi+1 = xi − ′
f (xi )
3x2i
starting with an initial guess at the root of x0 . This method only uses the four operations
given. Rather than do this “by hand” I’ll just use the python code newtonRaphson in the
script c 4 1 p 1.py. When we run that we get 4.21716332651 for the root estimate.
Problem 2
We first plot the function f (x) over a range of x values to get a feel for what this function
looks like. When we do that we get the plot given in Figure 1, from which we see three
real roots. The smallest positive real root seems to lie in the interval (0, 2). Using that as
a bracket of the root we can use the bisection method to refine this interval. We do this in
the python code c 4 1 p 2.py. When that code is run we get the value of 1.22999999998
for the root estimate.
22
2
1
0
f(x)
−1
−2
−3
−4
−5
−2.0
−1.5
−1.0
−0.5
0.0
x
0.5
1.0
1.5
2.0
Figure 2: A plot of the function f (x) for Problem 6.
Problem 3
I didn’t find the python module brent.py included in the source code from the book. Since
this module is given in the text it was easy to copy and paste that routine into a file and
then set up a calling routine to call this method. See the python code c 4 1 p 3.py where
I do this. When that code is run we get the value of 4.73004074486 for the root estimate.
Problem 4
See the python code c 4 1 p 4.py. When that code is run we also get the value of 4.73004074486
for the root estimate.
Problem 5
See the python code c 4 1 p 5.py. When that code is run we get the value of 7.068359375
for the root estimate.
Problem 6
We first plot our function over a range of x values to get a feel for were the roots are located.
When we do that we get the plot given in Figure 2. From that plot it looks like the leftmost
root is located in side (−1, 0) and the rightmost root is located inside (1, 2). Using these
23
500
400
300
f(x)
200
100
0
−100
−200
−300
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
7.5
8.0
x
Figure 3: A plot of the function f (x) for Problem 8. Note that for x < 4 the function f (x) is
very flat in that region. The tangent line through (x0 , f (x0 )) for x0 = 4 is drawn in red. The
next approximate root would be taken to be where this line intercepts the x-axis. Notice
that this intersection is “far to the right” of the root we seek and further iterations have
little chance of converging to the desired root.
intervals we can use the newtonRaphson routine to refine these estimates. See the python
code c 4 1 p 6.py. When that code is run we get the following
Leftmost root= -0.564326569397
Rightmost root= 1.20782767819
Problem 7 (the secant formula)
The secant method is implemented in the python code c 4 1 p 7.py. The implementation
given here was modeled after newtonRaphson in that it ensures the root is bracketed at each
iteration. When we run that code we get
Leftmost root= -0.56432657481
Rightmost root= 1.20782767948
Problem 8
Part (a-b): See the python code c 4 1 p 8.py. When that code is run we get the plot of
our function f (x) in Figure 3. As claimed, the smallest positive zero does lie in the interval
(4, 5) (another zero can be found in the interval (7.5, 8.0)).
24
1.8
1.6
1.4
1.2
f(x)
1.0
0.8
0.6
0.4
0.2
0.0
1.4
1.6
1.8
2.0
x
2.2
2.4
2.6
Figure 4: A plot of the function f (x) for Problem 9.
Problem 9
This problem is implemented in the the python code c 4 1 p 9.py. One small difficulty
with this problem is that since the root x = 2 is located at the bottom of a parabolic bowl
(see Figure 4) we have to have an initial bracket around the root that is quite small. This
is exactly what the routine rootsearch is for (to find a bracketing interval). When we run
the given python code we get the following
a, b (from rootsearch)=
root= 2.1
f(root)= 0.0
2.0 2.1
Note that if you know your root is a double root before you start it would be better use the
“multiple root” version of the Newton-Raphson method i.e.
xi+1 = xi − m
f (xi )
.
f ′ (xi )
where m is the roots multiplicity (two in this case). The text has a small comment about
this.
Problem 10
See the python code c 4 1 p 10.py. When that code is run we get the plot of our function
f (x) in Figure 5 and estimates of the three roots (using the brent solver) given by
25
8
6
4
2
f(x)
0
−2
−4
−6
−8
−10
−6
−4
−2
0
x
2
4
6
Figure 5: A plot of the function f (x) for Problem 10 over the interval (−6, +6).
The roots are:
-4.71238898039
-3.20883873198
1.57079632679
Problem 11
For this problem also see the python code c 4 1 p 10.py where when we run it using the
Newton-Raphson solver we get the following root estimates
-4.71238898041
-3.20883873198
1.5707963268
which are the same roots found in Problem 10.
Problem 12
For this problem also see the python code c 4 1 p 12.py where when we run it we get the
plot given in Figure 6. There we see that there are two roots: one between (−4, −2) and
one between (1, 3). Using the Newton-Raphson method we compute
Leftmost root=
-3.0
26
700
600
500
f(x)
400
300
200
100
0
−100
−6
−4
−2
0
x
2
4
6
Figure 6: A plot of the function f (x) for Problem 12.
Rightmost root=
2.1
Problem 13
For this problem also see the python code c 4 1 p 13.py where when we run it we get the
plot given in Figure 7. There we see that there are two positive roots: one between (0.5, 1.0)
and one between (1.5, 2). Using the Newton-Raphson method in each of these regions we
compute
Leftmost root= 0.791287847478
Rightmost root= 1.61803398876
Problem 14
For this problem also see the python code c 4 1 p 14.py where when we run it we get the
plot given in Figure 8. There we see that there are three positive roots: one between (1.5, 5.0),
one between (6, 7.5), and one between (7.5, 10). Using the Newton-Raphson method in each
of these regions we compute
First root= 2.85234189445
Second root= 7.0681743581
Third root= 8.42320393236
27
80
70
60
50
f(x)
40
30
20
10
0
−10
−1.0
−0.5
0.0
0.5
1.0
x
1.5
2.0
2.5
3.0
Figure 7: A plot of the function f (x) for Problem 13.
1.0
0.5
0.0
f(x)
−0.5
−1.0
−1.5
−2.0
−2.5
−3.0
0
5
10
x
15
20
Figure 8: A plot of the function f (x) for Problem 14.
28
200
150
f(x)
100
50
0
−50
1
2
4
3
5
6
x
Figure 9: A plot of the function f (x) for Problem 15.
Problem 15
For this problem see the python code c 4 1 p 15.py where when we run this code we get
the plot given in Figure 9 and the output
f_1=
f_2=
2.51660711564
15.7713075816
for the lowest two frequencies in units of Hz. Plotting our function is very helpful for it allows
us to visually see where the two roots are located. In that case we can then use the routine
newtonRaphson with brackets that allow us to specify each root of interest. Note that in
1
computing these we need to compute the moment of inertial for this beam as I = 12
wh3 .
Problem 16
For this problem see the python code c 4 1 p 16.py where when we run this code we get
the plot given in Figure 10 and the output
beta= 0.763400797561
cosh(beta)= 1.3058195668
for the solution for β and the value of cosh(β). Plotting our function is very helpful for
it allows us to visually see where the root in β is located. We can then use the routine
29
0.8
0.7
0.6
0.5
f(x)
0.4
0.3
0.2
0.1
0.0
−0.1
−0.5
0.0
0.5
1.0
1.5
2.0
x
Figure 10: A plot of the function f (x) for Problem 16.
newtonRaphson with brackets that allow us to specify the root of interest. To obtain the
value of σmax we would need to multiply cosh(β) by σ0 (which was not given a numerical
value in the problem).
Problem 17
For the given numbers in the python code c 4 1 p 17.py and for the given value σmax =
120 × 106 we plot the secant formula function
√ f (˜
σ) = σ˜ 1 + C1 sec C2 σ
˜ − σmax ,
as a function of σ˜ with C1 and C2 are constants determined by the formula given in the
book. When we do that we get Figure 11 and the output
C1, C2= 0.716623685777 9.38233281301e-05
sigma_tilde= 61083430.6544
for the solution for σ
˜ . Given this value of σ
˜ we would multiply by A = 25800 to determine
the value of P .
Problem 18
See the python code c 4 1 p 18.py for an implementation of this problem. When we run
that code we get Figure 12. In that figure we see two roots or possible values of the variable
30
3.0 1e8
2.5
2.0
1.5
f(x)
1.0
0.5
0.0
−0.5
−1.0
−1.5
0.0
0.2
0.4
0.6
0.8
x
1.0
1.4
1.2
1.6
1e8
Figure 11: A plot of the function f (x) for Problem 17.
5
4
f(x)
3
2
1
0
−1
0.0
0.2
0.4
0.6
0.8
1.0
x
Figure 12: A plot of the function f (x) for Problem 18.
31
12000
10000
8000
f(x)
6000
4000
2000
0
−2000
0.00
0.05
0.10
0.15
0.20
0.25
x
Figure 13: A plot of the function f (x) for Problem 19.
h. Running the newtonRaphson code with brackets specified for the output of the python
code
LHS, C1= 0.662923950114 0.022652622041
h_left= 0.26475526339
h_right= 0.49575512424
Problem 19
See the python code c 4 1 p 19.py for an implementation of this problem. When we run
that code we get Figure 13 and the output
t=
0.026322719651
Problem 20
See the python code c 4 1 p 20.py for an implementation of this problem. When we run
that code we get Figure 14 and the output
T2/T1=
5.41254824141
32
0.3
0.2
f(x)
0.1
0.0
−0.1
−0.2
−0.3
0.0
0.2
0.4
0.6
0.8
1.0
x
Figure 14: A plot of the function f (x) for Problem 20.
50
0
−50
f(x)
−100
−150
−200
−250
−300
0
2
4
6
8
10
12
14
x
Figure 15: A plot of the function f (x) for Problem 21.
33
10000
8000
f(x)
6000
4000
2000
0
−2000
0.55
0.60
0.65
0.70
0.75
x
0.80
0.85
0.90
0.95
Figure 16: A plot of the function f (x) for Problem 22.
Problem 21
I had to assume that there was a typo in this problem in expression for G should have been
G = −10 Joules otherwise normal ranges of T would not allow the logarithmic function to
obtain a value of G as large as specified in the book. With that condition see the python
code c 4 1 p 21.py for an implementation of this problem. When we run that code we get
Figure 15 and the output
T=
4.90242037274
Problem 22
See the python code c 4 1 p 22.py for an implementation of this problem. When we run
that code we get Figure 16 and the output
xi=
0.81712188891
Problem 23
See the python code c 4 1 p 23.py for an implementation of this problem. When we run
that code we get Figure 17 and the output
Leftmost intersection= [ 0.27942331
1.01961554]
34
5
4
3
2
1
0
−1
−2
−2
0
x
2
4
Figure 17: A plot of the two circles specified in Problem 23.
Rightmost intersection= [ 1.72057669
1.98038446]
Problem 24
This problem is implemented in the python code c 4 1 p 24.py. When we run that code
we get the following
[ 1.20782768
0.58842431]
for the solution near the point (1, 1).
Problem 25
See the python code c 4 1 p 25.py for an implementation of this problem. When we run
that code we get Figure 18 and the output
intersection 1=
intersection 2=
[ 0.31223443
[ 1.41423939
0.32279297]
6.33518146]
Notice that there are two solutions corresponding to y increased by 2π.
35
14
12
tan
arcsin(k=0)
arcsin(k=1)
10
8
6
4
2
0
0.0
0.2
0.4
0.6
0.8
x
1.0
1.2
1.4
1.6
Figure 18: A plot of the two circles specified in Problem 25.
Problem 26
Based on the problem statement the function f we seek the roots of is given by


(8.21 − a)2 + (0.0 − b)2 − R2
f (a, b, R) =  (0.34 − a)2 + (6.62 − b)2 − R2  .
(5.96 − a)2 + (−1.12 − b)2 − R2
See the python code c 4 1 p 26.py for an implementation of this problem. When we run
that code we get the output
a0, b0, R0= [4.8366666666666669, 1.8333333333333333, 4.5222000492029402]
a, b, R= [ 4.83010565 3.96992168 5.21382431]
Which correspond to our initial guess at the solution for a, b, and R and then the final
estimated values.
Problem 27 (an orbiting satellite)
Based on the problem statement, the function we seek to find the roots of is given by


C
− 6870
1+e sin(θ1 +α)

C
− 6728 
f (C, e, α) =  1+e sin(θ
.
2 +α)
C
−
6615
1+e sin(θ3 +α)
36
90°
135°
45°
8000
7000
6000
5000
4000
3000
2000
1000
0°
180°
225°
315°
270°
Figure 19: A plot of the orbit found for Problem 27.
Here θi are the three angle measurements measured in radians (not degrees). In the python
code c 4 1 p 27.py we use newtonRaphson2 to find estimate of C, e, and α. When we run
that code we get the output
initial values of C,e,alpha= [6737.666666666667, 0.00023933381867369263, 0.0]
final values of C,e,alpha= [ 6.81929379e+03 -4.05989591e-02
6.63142297e+01]
If we plot R as a function of θ we get Figure 19. Once we have values for the parameters C,
e, and α we look for the value of θ that makes R(θ) a minimum. To find this value of θ we
need to solve
C
R′ (θ) = −
(e cos(θ + α)) = 0 ,
(1 + e sin(θ + α))2
for θ. This means that
cos(θ + α) = 0 .
There are two solutions given by
π
− α mod 2π = 250.474515192 degrees
θ1 =
2
3π
θ2 =
− α mod 2π = 70.4745151918 degrees .
2
From Figure 19 we see that θ1 corresponds to the point with the largest radius and θ2 the
point with the smallest radius. See the python code c 4 1 p 27.py for an implementation
of this problem.
37
Problem 28 (a projectile)
From what we are given we know that the velocity of our projectile has components
x˙ = v cos(θ)
y˙ = −gt + v sin(θ) .
At the end of the trajectory to have the “angle of attack” specified we want to have
y˙
= −1 .
tan
x˙
(3)
Two other constraints specify the (x, y) location of the projectile at the end of the trajectory.
These equations are
x(t) = v cos(θ)t = 300
1
y(t) = − gt2 + v sin(θ)t = 61 .
2
For an initial guess at the solution we first choose θ0 =
in Equations 4 and 4 we need to have
π
4
(4)
(5)
which means that using this value
vt
√ = 300
2
vt
1 2
− gt + √ = 61 .
2
2
We can put the first equation in the second and solve for t. Using this we then get an initial
value for v. This problem is worked in the python code c 4 1 p 28.py. When we run that
code we get
initial values of v,theta,t= [60.779455899294014, 0.7853981633974483, 6.9803860932038475]
final values of v,theta,t= [ 60.02608349
0.8727932
7.7764307 ]
Problem 29
We need to measure θ3 in radians where 75 degrees is 1.308997 radians. The function we
want to find the zeros of has two components given by
150 cos(θ1 ) + 180 cos(θ2 ) − 200 cos(θ3 ) − 200 = 0
150 sin(θ1 ) + 180 sin(θ2 ) − 200 sin(θ3 ) = 0 .
Visually from the diagram it looks like the two roots would correspond
• θ1 ≈ θ3 and θ2 “small” say less than 45 degrees.
• θ1 “small” and θ2 “large” say greater than 45 degrees.
38
With initial guesses that mirror this sentiment, this problem is worked in the python code
c 4 1 p 29.py. When we run that code we get
initial values of theta_1,theta_2 (degrees)=
final values of theta_1,theta_2 (degrees)= [
initial values of theta_1,theta_2 (degrees)=
final values of theta_1,theta_2 (degrees)= [
[ 75. 10.]
54.98121678 23.00310093]
[ 10.
67.5]
20.01878322 51.99689907]
Problem Set 4.2
Notes on deflation of polynomials
We start with a root r of the nth degree polynomial Pn (x) and seek to write it in terms of
a n − 1 degree polynomial Pn−1 (x) such that Pn (x) = (x − r)Pn−1 (x). The left-hand-side of
this expression is
Pn (x) = a0 + a1 x + a2 x2 + · · · + an−1 xn−1 + an xn ,
while the right-hand-side of this expression is
(x − r)Pn−1 (x) = (x − r)(b0 + b1 x + b2 x2 + · · · + bn−1 xn−1 )
= b0 x + b1 x2 + · · · bn−2 xn−1 + bn−1 xn
− rb0 − rb1 x − · · · − rbn−2 xn−2 − rbn−1 xn−1
= −rb0 + (b0 − rb1 )x + (b1 − rb2 )x2 + · · ·
+ (bn−3 − rbn−2 )xn−2 + (bn−2 − rbn−1 )xn−1 + bn−1 xn .
Setting the coefficients of x equal on each side of the above we get
bn−1 = an
bn−2 − rbn−1 = an−1
bn−3 − rbn−2 = an−2
..
.
b1 − rb2 = a2
b0 − rb1 = a1
⇒
⇒
⇒
⇒
bn−2 = an−1 + rbn−1
bn−3 = an−2 + rbn−2
b1 = a2 + rb2
b0 = a1 + rb1 .
Which shows the general relationship for polynomial deflation
bn−1 = an
and bi = ai+1 + rbi+1
for i = n − 2, n − 3, . . . , 1, 0
(6)
Problems 1-9
All of these problems can be worked by “hand” or with calls to the deflPoly function
in the polyRoots module. See the python code c 4 2 p 1 9.py where we verify our hand
calculations using deflPoly.
39
Problem 1
We have
b2 = 3
b1 = 7 + (−5)(3) = −8
b0 = −36 + (−5)(−8) = 4 .
Thus
P3 (x) = (x + 5)(4 − 8x + 3x2 ) .
We can check this by multiplying the factored polynomial together as
(x + 5)(4 − 8x + 3x2 ) = 4x − 8x2 + 3x3 + 20 − 40x + 15x2 = 20 − 36x + 7x2 + 3x3 ,
which is P3 (x) as it should be.
Problem 2
We have
b3
b2
b1
b0
=1
= 0 + 1(1) = 1
= −3 + 1(1) = −2
= 3 + 1(−2) = 1 .
Thus
P4 (x) = (x − 1)(x3 + x2 − 2x + 1)
Problem 3
We have
b4
b3
b2
b1
b0
=1
= −30 + 6(1) = −24
= 361 + 6(−24) = 217
= −2178 + 6(217) = −876
= 6588 + 6(−876) = 1332 .
40
Problem 4
We have
b3
b2
b1
b0
=1
= −5 + 2i(1) = −5 + 2i
= −2 + (2i)(−5 + 2i) = −6 − 10i
= −20 + (2i)(−6 − 10i) = −12i .
Problem 5
We have
b2 = 3
b1 = −19 + ((3 − 2i)(3) = −10 − 6i
b0 = 45 + (3 − 2i)(−10 − 6i) = 3 + 2i .
Problem 6
We have
b2 = 1
b1 = 1.8 − 3.3(1) = −1.5
b0 = −9.01 + (−3.3)(−1.5) = −4.06 .
Problem 7
We have
b2 = 1
b1 = −6.64 + 0.64(1) = −6
b0 = 16.84 + 0.64(−6) = 13 .
Problem 8
We have
b2 = 2
b1 = −13 + (3 − 2i)(2) = −7 − 4i
b0 = 32 + (3 − 2i)(−7 − 4i) = 3 + 2i .
41
Problem 9
We have
b3
b2
b1
b0
=1
= −3 + (1 + 3i)(1) = −2 + 3i
= 10 + (1 + 3i)(−2 + 3i) = −1 − 3i
= −6 + (1 + 3i)(−1 − 3i) = 2 − 6i .
Problems 10-15
All of these problems can be worked with calls to the polyRoots function. See the python
code c 4 2 p 10 15.py. When we run that code we get (the rather ugly)
P 10:
[ 9.54249091e-18-1.j
1.10000000e+00+0.j
0.00000000e+00+1.j
-3.20000000e+00+0.j]
P 11:
[
1.+0.j
-1.+0.j
2.+0.j
-2.+0.j 156.+0.j]
P 12:
[ 1.+0.j 2.-1.j -3.+0.j -3.-1.j 2.+1.j -3.+1.j]
P 13:
[ 0.50000000 +0.00000000e+00j 1.00000000 -5.00000000e-01j
-1.99999999 -5.00865413e-08j 1.00000000 +5.00000000e-01j
-1.00000000 +2.00000000e+00j -1.00000000 -2.00000000e+00j
-2.00000001 +5.00865411e-08j]
P 14:
[ 2.+0.j -1.+0.j 3.-1.j -2.+0.j 3.+1.j -1.-2.j -1.+2.j 4.+0.j]
P 15:
[ 2.00000000e+00+0.j
2.19634966e-15+2.j -4.71844785e-16-3.j
-7.00000000e+00+0.j]
Problem 16
See the python code c 4 2 p 16.py. When we run that code we get
[ -0.62301963-24.03024141j -0.62301963+24.03024141j
-11.37698037+61.35447281j -11.37698037-61.35447281j]
42
Chapter 5 (Numerical Differentiation)
Notes on the Text
Notes on Example 5.3
We drop perpendiculars from the points B and C to the line segment AD. Call the intersection of the perpendicular dropped from B and the line segment AD the point B ′ and the
intersection of the perpendicular dropped from C and the line segment AD the point C ′ .
Then the right triangle CC ′ D has legs of lengths
CC ′ = a cos(α) + b cos(β)
C ′ D = d − a sin(α) − b sin(β) .
Since CC ′ D is a right triangle we can use the Pythagorean theorem to express the length of
the hypotenuse c in terms of the lengths of the triangles two sides as
2
2
CC ′ + C ′ D = c2 .
Replacing what we know about the lengths of the segments CC ′ and C ′ D we have the above
is equal to
(a cos(α) + b cos(β))2 + (d − a sin(α) − b sin(β))2 = c2 ,
which is the expression given in the book.
Problems
Problem 1
We start with the Taylor series for f (x + h2 ) and f (x − h1 ) as
h22 ′′
h32 ′′′
f (x + h2 ) = f (x) + h2 f (x) + f (x) + f (x) + O(h42 )
2
3!
2
3
h
h
f (x − h1 ) = f (x) − h1 f ′ (x) + 1 f ′′ (x) − 1 f ′′′ (x) + O(h41 ) .
2
3!
′
If we multiply the first equation by h1 and the second equation by h2 and then add the two
we get the following
h1 f (x + h2 ) + h2 f (x − h1 ) = (h1 + h2 )f (x) +
1
h1 h22 + h2 h21 f ′′ (x)
2
1
+ (h1 h32 − h2 h31 )f ′′′ (x) + O(h1 h42 ) + O(h2 h41 ) .
6
43
Solving for f ′′ (x) in the above we get
2
[h1 (f (x + h2 ) − f (x)) + h2 (f (x − h1 ) − f (x))]
h1 h2 (h1 + h2 )
h31
h32
1
′′′
+O
.
− (h2 − h1 )f (x) + O
3
h1 + h2
h1 + h2
f ′′ (x) ==
If h1 ≈ h2 ≈ h then our truncation error is O(h2 ).
Problem 2
To begin recall that the first backwards finite difference approximation for f ′ (x) is given by
f ′ (x) =
1
(f (x) − f (x − h)) + O(h) ,
h
and the first backwards finite difference approximation for f ′′ (x) is given by
f ′′ (x) =
1
(f (x) − 2f (x − h) + f (x − 2h)) + O(h) .
h2
Now using the fact that f ′′′ (x) = (f ′′ (x))′ we have
d 1
′′′
f (x) =
(f (x) − 2f (x − h) + f (x − 2h))
dx h2
1
1 1
(f (x) − 2f (x − h) + f (x − 2h)) − 2 (f (x − h) − 2f (x − 2h) + f (x − 3h))
=
h h2
h
1
= 3 (f (x) − 3f (x − h) + 3f (x − 2h) − f (x − 3h)) .
h
Notice that the coefficients in this expression match the ones in the third row in Table 5.2b.
Problem 3
The central difference approximation for f ′′ (x) is given by
f ′′ (x) =
1
(f (x + h) − 2f (x) + f (x − h)) + O(h2 ) .
2
h
(7)
To apply Richardson extrapolation we need to write our target G ≡ f ′′ (x) in the form of
G = f ′′ (x) = g(h) + chp ,
so based on Equation 7 we need to define the expression g(h) be the terms involving f (·)
on the right-hand-side of Equation 7 and to take p = 2. Next we evaluate g(·) for a fixed h
and then at h2 and construct a more accurate estimate of G = f ′′ (x) by combining g(h) and
g h2 . Namely as p = 2 we have a more accurate estimate given by
22 g h2 − g(h)
4g h2 − g(h)
G=
=
.
22 − 1
3
44
With the expression for g(h) above we find that G is given by
4
1
h
G= g
− g(h)
3
2
3
h
1 1
h
4 4
− 2f (x) + f x −
−
f x+
(f (x + h) − 2f (x) + f (x − h))
=
3 h2
2
2
3 h2
16
16
32
h
h
1
2
1
= 2f x +
+ 2f x −
− 2 f (x) − 2 f (x + h) + 2 f (x) − 2 f (x − h)
3h
2
3h
2
3h
3h
3h
3h
1
10
1
16
h
16
h
= − 2 f (x − h) + 2 f x −
− 2 f (x) + 2 f x +
− 2 f (x + h) .
3h
3h
2
h
3h
2
3h
This expression for f ′′ (x) will be accurate to O(h4 ).
Problem 4
We want to derive the second forward finite difference approximation for f ′′′ (x) from Taylor
series. To do this we need to consider the Taylor expansions of f (x+h), f (x+2h), f (x+3h),
and f (x + 4h) and use these four equations to solve for f ′ (x), f ′′ (x), and f ′′′ (x). The
solution for f ′′′ (x) is what we want. As this involves a great deal of algebra using a algebraic
manipulation program (like Mathematica) would be helpful.
Problem 5
We want to derive the first central difference approximation for f (4) (x). To do this we will
write down the Taylor expansions of f (x + h), f (x − h), f (x + 2h), f (x − 2h) and then solve
them for f ′ (x), f ′′ (x), f ′′′ (x), and f (4) (x). The expression for f (4) (x) will be the desired
result.
Problem 6
From the data given we must use forward differences with h = 0.01 to estimate f ′ (2.36) and
f ′′ (2.36) . We find
1
(−3f (2.36) + 4f (2.37) − f (2.38)) + O(0.012)
2(0.01)
= 0.424 + O(0.0001) .
f ′ (2.36) =
For f ′′ (x) we will use
1
(2f (2.36) − 5f (2.37) + 4f (2.38) − f (2.39)) + O(0.012)
0.012
= −0.2 + O(0.0001) .
f ′′ (2.36) =
45
Problem 7
Notice that the spacing between the points x is not uniform and thus we cannot use the
formulas for the approximate derivatives specified on a uniform grid. Instead we will fit
a quadratic polynomial P2 (x) = a0 + a1 x + a2 x2 to the three given points and return the
values of the first and second derivatives of this polynomial. Following the example in this
chapter we find that the coefficients are given by a0 = 1.02597000, a1 = −0.06783333, and
a2 = −0.11666667. Using these we had
f ′ (1) ≈ P2′ (1) = a1 + 2a2 (1) = −0.3011667
f ′′ (1) ≈ P2′′ (1) = 2a2 = −0.2333333 .
Problem 8
Note that these points are evenly spaced with five points spaced at the smaller spacing of
h2 = 0.08 and three points spaced at twice that spacing of h1 = 0.16 = 2h2 . To evaluate
f ′′ (1) we will first use the second order central difference approximation with h1 = 0.16
1
(f (1 − 0.16) − 2f (1) + f (1 + 0.16))
0.162
1
(f (0.84) − 2f (1) + f (1.16)) = 0.3687109 .
=
0.162
Next we use a second order finite difference approximation with h2 = 0.08 as
f ′′ (1) =
f ′′ (1) =
1
(f (0.92) − 2f (1) + f (1.08)) = 0.3682812 .
0.082
Using Richardson’s extrapolation with p = 2 (since our methods are second order) we get
f ′′ (1) =
4f ′′h=0.08 (1) − f ′′ h=0.16 (1)
= 0.368138 .
3
In this case this is accurate to O(h41 ) = O(0.164) = O(0.00065536).
Problem 9
To obtain f ′ (0.2) as accurately as possible with the given data I will use Richardson’s extrapolation with a larger spacing of h1 = 0.2 and a smaller spacing of h2 = h21 = 0.1. When
h1 = 0.2 we have as our estimate of f ′ (0.2) the following
f ′ (0.2) =
1
1
(f (x + h1 ) − f (x − h1 )) =
(f (0.4) − f (0)) = 0.6124525 ,
2h1
2(0.2)
which is an O(h21 ) = O(0.04) approximation. Next we use h2 = 0.1 and central differences
again to get
f ′ (0.2) =
1
1
(f (x + h2 ) − f (x − h2 )) =
(f (0.3) − f (0.1)) = 0.57284 ,
2h2
2(0.1)
46
which is an O(h22 ) = O(0.01) approximation. Using Richardson’s extrapolation our estimate
(since p = 2) is given by
f ′ (0.2) =
4f ′ h=0.1 (0.2) − f ′ h=0.2 (0.2)
= 0.5596358 .
3
This approximation is O(h41 ) = O(0.24) = O(0.0016).
Problem 10
Part (a): For this part we will use
f ′ (x) =
1
(f (x + h) − f (x)) ,
h
for h ranging from small to large.
Part (b): For this part we will use
f ′ (x) =
f (x + h) − f (x − h)
2h
for h ranging from small to large. These approximations are implemented in the python
code c 5 1 p 10.py when we run that code we get the following output
First forward estimates:
h
error
1.000000
0.440217
0.500000
0.204307
0.250000
0.096467
0.100000
0.037007
0.050000
0.018307
0.010000
0.003707
0.001000
0.006707
0.000100 -0.003293
0.000050
0.096707
0.000010
0.696707
0.000001
0.696707
First central estimates:
h
error
1.000000
0.110447
0.500000
0.028667
0.250000
0.007247
0.100000
0.001157
0.050000
0.000307
0.010000 -0.000293
0.001000
0.001707
47
0.000100
0.000050
0.000010
0.000001
-0.003293
-0.003293
0.196707
0.696707
In this problem when h is large the error in the derivative is dominated by the truncation
error in our forward (or central) difference approximation. When h is small the truncation
error is small but the round-off error (in doing the subtraction) is poor. This is why the
error in each case initially starts to decrease as we make h smaller but at some point the
round-off error gets large and our approximation starts to get worse.
Problem 11
With the given data we can fit a quadratic P2 (x) = a0 + a1 x + a2 x2 and we find that
a0 = 21.441963, a1 = −6.171548, and a2 = −4.038080. Using these we can approximate
f ′ (0) and f ′′ (0) as
f ′ (0) ≈ a1 = −6.171548
f ′′ (0) ≈ 2a2 = −8.076161 .
Note that this polynomial does not fit the data very well so other methods may work better
(i.e. spline approximation techniques may work better).
Problem 12
We will evaluate x(θ) at the values of θ given by 0, 5, · · · , 175, 180 (converted to radians) to
get x(θi ). Since we want to approximate the second time derivative of x(t) note that
dx
dx dθ
dx
.
=
= 5000
dt
dθ dt
dθ
From this we have that the second derivative of x with respect to time t is given by
2
d2 x
2d x
=
5000
.
dt2
dθ2
2
We can evaluate the θ derivatives i.e. ddθx2 using forward/central/backwards differences. This
is implemented in the python code c 5 1 p 12.py. When we run that code we get the
2
following output representing ddθx2 at each of the 37 values of θ.
[-127.49886998 -125.0390675 -122.57926502 -118.51509775 -112.90079335
-105.81412576 -97.35804442 -87.66240481 -76.8854657
-65.21469813
-52.86633952 -40.08305933 -27.12912902 -14.2826618
-1.8248393
9.97343454
20.86639466
30.64816628
39.16545016
46.32624656
48
52.10323557
63.09519533
58.42762596
54.16221019
56.53130487
62.76202224
57.21777325
53.75075741]
59.69972693
61.9837691
56.13656311
61.74036124
60.91504984
55.2412398
62.81376971
59.69141392
54.57366298
These numbers would need to be multiplied by 50002 to convert them into second derivatives
of x with respect to time t.
Problem 13
We are given expressions for x any y in terms of the angles α and β which are themselves
functions of time t. We are asked to calculate the planes speed v and the climb angle γ
which are given in terms of the time derivatives of x and y as
p
v = x˙ 2 + y˙ 2
y˙
−1
γ = tan
.
x˙
To use these we need to relate the time derivative of x and y to the time derivatives of α
and β. We have
!
1
sec2 (β)β˙
1 dx
d
tan(β)
tan(β)(sec2 (β)β˙ − sec2 (α)α)
˙
=
x˙ =
=
−
a
a dt
dt tan(β) − tan(α)
tan(β) − tan(α)
(tan(β) − tan(α))2
− tan(α) sec2 (β)β˙ + tan(β) sec2 (α)α˙
(tan(β) − tan(α))2
tan(α) tan(β)
1 dy
d
1
y˙ =
=
a
a dt
dt tan(β) − tan(α)
2
sec (α) tan(β)α˙ + tan(α) sec2 (β)β˙
tan(α) tan(β)(sec2 (β)β˙ − sec2 (α)α)
˙
=
−
2
tan(β) − tan(α)
(tan(β) − tan(α)
sec2 (α) tan2 (β)α˙ − sec2 (β) tan2 (α)β˙
.
=
(tan(β) − tan(α)2
=
Next we estimate α˙ and β˙ using central differences and the data given and then use the
above formulas to estimate x˙ and y.
˙ We do this in the python code c 5 1 p 13.py. When
we run that code we get estimate of x˙ and y˙ at all three times given by
x_dot=
y_dot=
[-9.08948709 -8.56152742 -8.03011676]
[ 13.90177307 13.08554517 12.25522978]
One could then compute v and γ using the above formulas. Something must be incorrect
with the above calculations since the sign of x˙ should be positive. If anyone sees anything
incorrect that I have done please contact me.
49
Chapter 6 (Numerical Integration)
Problem Set 6.1
Problem 1 (the recursive trapezoidal rule)
In the python code c 6 1 p 1.py we implement the recursive trapezoidal rule using the
trapezoid function provided by the book. When we run that code we get
Integral= 0.272198261288
nPanels= 2
Indicating that we only needed two panels (three points) to integrate this integral.
Problem 2
We have P = F v and F = m dv
so P = mv dv
. Solving this for dt gives
dt
dt
mv
dv .
dt =
P
If we integrate both sides of this from t1 to t2 we get the expression quoted in the book
Z t2
v
∆t = m
dv .
t1 P
From the problem we are given v and P evaluated at points in time between zero and
six seconds and therefore we can use the trapezoidal rule in each panel to perform the
integration even though the spacing in ∆v is not uniform. The only wrinkly is that we will
have to approximate the total integral I using trapezoids of different sizes in each panel i.e.
I =m
7
X
Ii .
i=1
Where we will evaluate the integral over panel i or Ii with the trapezoidal rule as
(vi − vi−1 ) vi
vi−1
for 1 ≤ i ≤ 7 .
Ii =
+
2
Pi Pi−1
Note in the first interval I1 when we evaluate the value of Pv00 we get 00 an indeterminate
expression. To simplify things I’ll assume this evaluates to zero. With the number from the
text the first two Ii are given by
(1 − 0) 1.0
+0
I1 =
2
4.7
1.8 − 1 1.8
1.0
I2 =
.
+
2
12.2 4.7
50
In the python code c 6 1 p 2.py we finish the above calculations. When we run that code
we get the following
I=
0.755630597921
for the value of the integral (without the multiplication of the mass m).
Problem 3
See the python code c 6 1 p 3.py for an implementation of this procedure. When we run
that code we get the following output
nPanels= [2, 4, 6]
Integrals= [-0.66666666666666663, -0.66666666666666652, -0.66666666666666674]
Problem 4
To transform this integral to have fixed integration limits we let x3 = 1t which means that
x = t−1/3 . For this expression we have that dx = − 31 t−4/3 dt. The limits of the integral in x
transform to limits in t of
x = 1 becomes t = 1
x = +∞ becomes t = 0 .
Using these results our integral now becomes
Z
Z 1
1 1
1 −4/3
−4/3 −1
t
dt =
(1 + t4/3 )−1 dt .
I=
(1 + t
)
3
3
0
0
Note that this last integral has fixed integration limits which we can then integrate numerically with the trapezoidal rule. In the python code c 6 1 p 4.py we finish the above
calculations. When we run that code we get the following
Integral=
0.243698304044
Problem 5
I decided to integrate the given numerical values of x and F (x) function using the trapezoid
rule. See the python code c 6 1 p 5.py. When we run that code e get the following output
51
I= 74.4
v = 44.5421149026
The first number is the approximate integral of the force while the second number is the
velocity (in meters per second).
Problem 6
See the python code c 6 1 p 6.py, which is really just a call to the function romberg. When
we run that code e get the following output
(18.666666666666668, 8)
Problem 7
Note we can use Romberg integration (by hand) since we can compute Ri,1 (for 1 ≤ i ≤ 3)
using the composite trapezoid rule with n = 2i−1 panels. We can then combine these results
(and subsequent results) using Richardson extrapolation. Note that the exact formula used
for Richardson extrapolation changes as we move from column to column since the order
of the approximation increases. We follow this procedure in the python code c 6 1 p 7.py.
When we run that code we get the following output
[[ 3.14159265
[ 1.96349541
[ 1.52068792
0.
1.57079633
1.37308543
0.
]
0.
]
1.3599047 ]]
Problem 8
For the integral given
I=
if we let v =
√
Z
1
0
sin(x)
√ dx ,
x
√ and I can be written as
x (so that x = v ) then we have dv = 2dx
x
Z 1
I =2
sin(v 2 )dv .
2
0
This later integral has no indeterminacy at x = 0. To integrate this later expression see the
python code c 6 1 p 8.py. When we run that code we get the following output
(0.62053660344524453, 32)
52
Problem 9
From Section 3.3 in the book when modeling our function with cubic splines the approximate
function for all x such that xi ≤ x ≤ xi+1 is given by equation 3.10 where
ki+1 (x − xi )3
ki (x − xi+1 )3
− (x − xi+1 )(xi − xi+1 ) −
− (x − xi+1 )(xi − xi+1 )
fi,i+1 (x) =
6 xi − xi+1
6 xi − xi+1
yi (x − xi+1 ) − yi+1 (x − xi )
+
,
xi − xi+1
and i = 0, 1, 2, . . . , n − 1. The values of the curvatures, ki are given by the solution of a
tridiagonal system that depends on the pairs (xi , yi ) . When xi+1 − xi = h i.e. uniform
spacing we have that
ki
(x − xi+1 )3
(x − xi )3
ki+1
fi,i+1 (x) =
−
−
+ h(x − xi+1 ) −
+ h(x − xi+1 )
6
h
6
h
−yi (x − xi+1 ) + yi+1 (x − xi )
+
.
h
If we integrate this expression over [xi , xi+1 ] we get
Z
xi+1
fi,i+1 (x)dx =
xi
−
=
=
x
x
ki
(x − xi+1 )4 h(x − xi+1 )2 i+1 ki+1
(x − xi )4 h(x − xi+1 )2 i+1
−
+
−
−
−
6
4h
2
6
4h
2
xi
xi
xi+1
xi+1
2
2
yi (x − xi+1 ) yi+1 (x − xi ) +
2h
2h
xi
xi
4
4
h
3
h
h3
ki+1
h
yi 2 i yi+1 2
ki h
−
−0+
h +
h
−
−
− 0−
6 4h
2
6
4h
2
2h
2h
h3
h
−(ki + ki+1 ) + (yi + yi+1 ) ,
24
2
when we simplify. Given this expression we can sum for i = 0, 1, 2, . . . , n − 1 to get an
approximation to our integral I. Doing this we have
n−1
n−1
hX
h3 X
I≈
(yi + yi+1 ) −
(ki + ki+1 )
2 i=0
24 i=0
=
h
h3
(y0 + 2y1 + 2y2 + · · · + 2yn−1 + yn ) − (k0 + 2k1 + 2k2 + · · · + 2kn−1 + kn ) ,
2
24
which is the desired result.
Problem 10
For the integral given
I=
Z
0
π/4
dx
p
,
sin(x)
53
if we let sin(x) = t2 then taking the differential of both sides we get
cos(x)dx = 2tdt .
Thus solving for dx we get
dx =
2tdt
2tdt
2tdt
=√
=p
.
2
cos(x)
1 − t4
1 − sin (x)
Using this our integral becomes
Z ( 1 )1/4
Z ( 1 )1/4 2
2
1
dt
2tdt
√
√
.
=2
I=
t
1 − t4
1 − t4
0
0
This integral does not have a singularity at 0. To integrate this later expression see the
python code c 6 1 p 10.py. When we run that code we get the following output
(1.7911613389539645, 64)
Problem 11
See the python file c 6 1 p 11.py where this problem is worked. When we run that function
we get the following
pi/2= 1.57079632679
theta0 values= [0, 15, 30, 45]
h(theta0)= [1.5707963267948966, 1.5775516530701426, 1.598142002573657, 1.6335863090850953]
Problem 12
See the python file c 6 1 p 12.py where this problem is worked and we use the romberg
function to integrate this function. When that code is run we get the following output from
this function
(0.40629888639968503, 32)
Problem 13
The expression for the singularity (1 − z 2 )−1/2 is the kernel in Gauss-Chebyshev integration
R +1
but I don’t see how to transform the integral into the required form −1 (1 − x2 )−1/2 f (x)dx
for that technique. Instead since Gauss-Legendre integration does not evaluate the integrand
at the limits of integration we will use that technique. See the python file c 6 1 p 13.py
where this problem is worked and we use the romberg function to integrate f (x) and then
compute v0 . When that code is run we get the following output from this function
54
2.49767483249
Problem 14
There are two difficulties with this problem as written. One is that the the point u = 0
should really be evaluated separately. In that case the upper limit of the integrand is +∞
and since the integral converges when we multiply by u = 0 we have that g(0) = 0. The
second is that the integrand has a removable singularity at x = 0 which gives the standard
methods introduced in this chapter trouble. We can find the value that the integrand has
when x = 0 by repeated application of L’ Hopital’s rule. We have
4x3 ex + x4 ex
x4 ex
=
lim
since our initial limit was of the type 0/0
x→0 2(ex − 1)ex
x→0 (ex − 1)2
4x3 + x4
canceling ex in top and bottom
= lim
x→0 2(ex − 1)
12x2 + 4x3
= lim
again our limit was of the type 0/0
x→0
2ex
= 0.
lim
See the python file c 6 1 p 14.py where this problem is worked. When we run this script
we get the following output
[ 0.
0.05
0.6
0.65
[ 0.
0.2025676
0.29126528
0.31363118
0.1
0.15
0.7
0.75
0.00324692
0.22885759
0.29699209
0.31557282
0.2
0.25
0.8
0.85
0.02527367
0.24861767
0.30165145
0.31724418]
0.3
0.35
0.9
0.95
0.07099747
0.26360757
0.30548685
0.4
0.45
1. ]
0.12287827
0.27513603
0.30867793
0.5
0.55
0.16768644
0.28413572
0.31135885
Here the top array represents the value of u and the bottom array represents the value of
g(u).
Problem 15
Since the integrand of this integral decreases exponentially we will evaluate the infinite
integral by truncating it at some finite value T as
Z ∞
Z T
Z ∞
2t
2t
2t
2 −2t/t0
2
2 −2t/t0
2
2 −2t/t0
2
i0 e
sin
dt =
dt +
dt .
i0 e
sin
i0 e
sin
t0
t0
t0
0
0
T
Then to evaluate the total integral to an accuracy of ǫ we need to pick T such that
Z ∞
2t
2 −2t/t0
2
≤ ǫ.
dt
i
e
sin
0
t0
T
55
We can bound the right-hand-side as
Z ∞
Z ∞
2t
2
2 −2t/t0
2
dt ≤ i0
e−2t/t0 dt
i0 e
sin
t
0
T
∞
T
−2t/t0 t0 −2T /t0
2e
2
= i0 0 + e
= i0
.
2
− t2 0
T
We thus want to pick a value of T such that
i20 t0 −2T /t0
e
≈ ǫ,
2
or
2ǫ
t0
.
T = − log 2
2
i0 t0
Our true integral would be multiplied by R which would reduce the value of ǫ by a factor
of R. With ǫ = 10−6 , i0 = 100, and t0 = 0.01 the above gives T = 0.0851719319142. Using
this upper limit, our integral is then performed in the python code c 6 1 p 15.py. When
we run that code we get the following
T= 0.0851719319142
Our integral estimated using the recursive trapezoid rule...
Integral= 9.99999930184
nPanels= 1024
Our integral estimated using romberg...
(9.9999993172368455, 256)
Problem Set 6.2 (Gaussian Integration)
Problem 1
Lets first check to see if there are any singularities in this integral. First note that ln(x) is
singular at x = 0 and complex for x < 0. This is not a problem here since the integration
region is (1, π). Next the denominator of the fraction is the polynomial x2 − 2x + 2 which
has zeros at
p
√
2 ± 4 − 4(2)
2 ± −4
=
= 1 ± i.
2
2
which are complex and again don’t represent a problem. For this we will use the routine
gaussQuad. See the python code c 6 2 p 1.py where this problem is worked. When we run
that code we get the following
Number
Number
Number
Number
of
of
of
of
nodes=
nodes=
nodes=
nodes=
2;
4;
6;
8;
I=
I=
I=
I=
0.606725
0.584768
0.584943
0.584943
56
Problem 2
Note that (1 − x2 )3 is a sixth order polynomial thus using Gauss-Laguerre integration we
can integrate this integral exactly when 2n + 1 ≥ 6 so n ≥ 2.5. If we take n = 3 and using
Table 6.4 from the book we can evaluate this integral. See the python code c 6 2 p 2.py
where this problem is worked. When we run that code we get the following
I (Gauss-Laguerre)= -652.933958
Problem 3
For the given integral
I=
Z
π/2
0
dx
p
,
sin(x)
if we let sin(x) = t2 with a differential of cos(x)dx = 2tdt. Solving for dx there we have
dx =
2tdt
2tdt
2tdt
√
=
=p
.
cos(x)
1 − t4
1 − sin2 (x)
Under this change of variable the limits of our integral transform as x = 0 ⇒ t = 0 and
x = π2 ⇒ t = 1. With these I becomes
1
Z 1
1 2t
1
√
p
I=
dt
dt = 2
4
2
1−t
(1 + t )(1 − t2 )
0 t
0
Z 1
Z 1
1
1
1
1
√
√
√
√
dt =
dt ,
=2
1 − t2
1 + t2
1 − t2
1 + t2
−1
0
Z
1
which is a Gauss-Chebyshev integral with the function f (t) = √1+t
2 integrated against the
2 −1/2
Chebyshev kernel (1 − t )
. Because of this to integrate we will use Gauss-Chebyshev
integration. See the python code c 6 2 p 3.py where this problem is worked. When we run
that code we get the following
Gauss-Chebyshev (6 nodes); I_approx=
2.622027; error=
Problem 4
Note that the exact value of this integral is given by
Z π
sin(x)dx = − cos(x)|π0 = −(1 − 1) = 2 .
0
57
0.000033
The truncation error in using Gauss-Legendre integration is given in Eq. 6.26 (or Eq. 6.30
for more general limits of integration). This later equation is given by
En =
(b − a)2n+3 [(n + 1)!]4 (2n+2)
f
(c) for a < c < b .
(2n + 3)[(2n + 2)!]3
(8)
To integrate with four notes we use n = 3 (since some of the points are symmetrical i.e. they
come in pairs). For this problem f (x) = sin(x) and thus all derivatives of this will be either
± sin(x) or ± cos(x) both of which are less than one in magnitude. If we take f (2n+2) (c) at
its upper bound of one we find
E3 ≤ 1.676447 10−5 .
Problem 5
We can evaluate the given integral in the following way (using integration by parts)
Z ∞
I=
e−x sin(x)dx
0
Z ∞
−x ∞
cos(x)e−x dx
= − cos(x)e 0 −
Z ∞ 0
= −(0 − 1) −
cos(x)e−x dx
0
Z ∞
∞
−x
−x
= 1 − e sin(x) 0 +
sin(x)e dx
0
= 1 − (0 − 0 + I) .
Solving for I we have that I = 21 for the exact value. For Gauss-Laguerre quadrature the
truncation error is given by Eq. 6.35 in the book or
En =
[(n + 1)!]2 (2n+2)
f
(c) .
(2n + 2)!
(9)
Here f (x) = sin(x) and thus all derivatives of this will be either ± sin(x) or ± cos(x) both of
which are less than one in magnitude. We will have the desired accuracy in our numerical
integral if we pick an n such that
|En | ≤
[(n + 1)!]2
≤ 10−6 .
(2n + 2)!
The value of n = 11 gives |En | ≤ 3.698012 10−7. Thus we need n + 1 = 12 nodes.
Problem 6
so that dx = dt2 and notice that the argument of the
To simplify this integral we let x = 1+t
2
square root is now
1+t
1+t
1
x(1 − x) =
1−
= (1 − t2 ) .
2
2
4
58
With this I becomes
Z +1
Z +1
Z 1
(1 + t + 1) dt
t+2
(2x + 1)
√
p
√
dx =
=
I=
dt ,
2
1−t
2
1 − t2
x(1 − x)
−1
−1
0
2
which is a Gauss-Chebyshev integral with the function f (t) = t + 2 (a first order polynomial)
integrated against the Chebyshev kernel (1 − t2 )−1/2 so we can evaluate this exactly with
n = 0 (i.e. one node) using Gauss-Chebyshev integration. Gauss-Chebyshev integration
(Eq. 6.1) is
Z
+1
−1
2 −1/2
(1 − x )
n
(2i + 1)π
π X
f (xi ) with xi = cos
.
f (x)dx ≈
n + 1 i=0
2n + 2
(10)
Taking n = 0, the integral we want (where f (t) = t + 2) becomes
Z +1
π π
= πf (0) = 2π ,
(1 − x2 )−1/2 f (x)dx ≈ f (x0 ) = πf cos
1
2
−1
which is actually exact as argued above.
Problem 7
This integral is singular at x = 0 due to the logarithmic factor. This is not a problem
for Gauss-Legendre quadrature since the integrand is never evaluated at the end points of
integration. Thus as a first attempt we will use Gauss-Legendre quadrature to solve this
problem. To figure out how many nodes we need to obtain the desired accuracy we could try
to make En in Equation 8 less than 10−4 . In this problem f (x) = sin(x) ln(x) and to bound
En we need to bound derivatives of f (x) for various powers. For the first derivative we find
f (1) (x) =
sin(x)
+ cos(x) ln(x) ≤ 1 + ln(π) ,
x
but we would need to find bounds on much higher derivatives to evaluate En in for larger n
to use Equation 8 and this is difficult/tedious. Another method to approximate this integral
would be to use Gauss quadrature with a logarithmic singularity for part for the part of the
integral. That is we could write our integral as
Z π
Z 1
Z π
sin(x) ln(x)dx =
sin(x) ln(x)dx +
sin(x) ln(x)dx .
0
0
1
The first integral can be integrated using Gaussian integration with a logarithmic singularity
and the second can be integrated using Gaussian-Legendre integration. What’s nice about
the first integral is that the function we are integrating the logarithmic kernel against is sin(x)
which has all derivatives bounded by one. Thus from the books Eq. 6.39 the truncation error
E1 (the 1 is for the first integral) using this technique is
E1 =
k(n)
f (2n+1) (c) ,
(2n + 1)!
59
(11)
which for the values of k(n) given in the book and the function f (x) = sin(x) here gives
|E1 | ≤
k(n)
1
<
.
(2n + 1)!
(2n + 1)!
We would then want to pick n for this first integral to guarantee that |E1 | ≤ 12 10−4 . For
the second integral in this splitting we would be back to using Equation 8 to bound the
truncation error for the in such a way that that |E2 | ≤ 21 10−4 . Unfortunately to bound
|E2 | we again have to be able to bound derivatives of f (x) = sin(x) ln(x) which could be
difficult/tedious.
In the end it might just be easier to use Gauss-Legendre integration on the full integral as
given but increasing the number of nodes used in its estimation it until we get the desired
relative accuracy of four decimal places. This later approach is worked in the python code
c 6 2 p 7.py. When we run that code we get the following
(0.64103922896858212, 11)
which is the integrals value and the number of nodes used in computing it.
Problem 8
We first derive the true value of the integral as follows
Z π
I=
x sin(x)dx
0
Z π
π
= − x cos(x)|0 +
cos(x)dx
0
= −(−π) + sin(x)|π0 = π .
Then the desired procedure is implemented in the python code c 6 2 p 8.py. When we run
that code we get the following
Gauss-Legendre (
3 nodes); I=
3.143774, true error=
-0.002182
Problem 9
This problem is implemented in the python code c 6 2 p 9.py. When we run that code we
get the following
(2.5015673863717733, 4)
which is the integrals value and the number of nodes used in computing it.
60
Problem 10
For our integral
I=
Z
0
∞
xdx
,
+1
ex
we let e = so that the differential of this expression is given by ex dx = − t12 dt. We can
write this as
1
dt
dx = t − 2 dt = − ,
t
t
to get
Z 0
Z 1
ln(t) − dtt
ln(t)
dt .
I=
− 1
=−
+1
1
0 1+t
t
x
1
t
To evaluate this we can use Gauss-Quadrature with a logarithmic singularity
Z 1
n
X
f (t) ln(t)dt ≈ −
Ai f (ti ) ,
0
i=0
with the values of Ai and ti given in Table 6.6 from the book. The function we are integrating
1
. By taking derivatives of f (t) we can conclude that
against ln is f (t) = − 1+t
f (n) (t) =
(−1)n−1 n!
.
(1 + t)n+1
Thus on the domain 0 < t < 1 we have that
|f (n) (t)| ≤ n! ,
is a bound for the derivative of the function f (t). This is not a very good bound since the
truncation error for Gauss quadrature with a logarithmic singularity has a truncation error
given by the books Eq. 6.39 or
E=
k(n)
f (2n+1) (c) for 0 < c < 1 ,
(2n + 1)!
(12)
if we use the bound derived above on our derivative of f this becomes for this particular
problem
k(n)
(2n + 1)! = k(n) .
|E| ≤
(2n + 1)!
The book only gives values of the function k(n) for n = 1, n = 2, and n = 3. When n = 3 we
have that |E| ≤ k(3) = 0.00001 = 10−5 which seems close enough to the desired accuracy of
10−6 that we should use it. It looks like k(n) is a decreasing sequence in n so one could use a
larger value of n i.e. n > 3 if one desired more accuracy. This problem is then implemented
in the python code c 6 2 p 10.py. When we run that code we get the following
I=
0.822466107543
which is the integrals approximate value.
61
Problem 11
We start with the integral
S=2
Z
a
−a
s
1+
dy
dx
2
dx .
Taking the derivative with respect to x of
x2 y 2
+ 2 = 1,
a2
b
we get
2x 2y dy
+ 2
= 0.
a2
b dx
If we solve for
dy
dx
we get
b2 x
dy
=− 2 .
dx
a y
We then square that expression to get
From
x2
a2
+
y2
b2
dy
dx
2
=
b4 x2
.
a4 y 2
= 1 we can solve for y 2 and put this in the above to get
dy
dx
2
b2
= 2
a
x2
a2 − x2
.
If we put this into the expression for S we get
Z as
b2
x2
1+ 2
S=2
dx .
a
a2 − x2
−a
Notice that this integral is singular at x = ±a. This is not a problem for Gauss-Legendre
quadrature since it does not actually evaluate the integrand at the limits of the integration.
We implement this integration in the python code c 6 2 p 11.py. When we run that code
we get the following output
I=
9.670165 with m=
190 nodes
Problem 12
This problem is implemented in the python code c 6 2 p 12.py. When we run that code
we get the following
62
0.1 (0.11246291602487697, 3)
0.25 (0.27632639016449356, 4)
0.5 (0.52049987679724563, 4)
1.0 (0.84270078611625499, 5)
2.0 (0.99532227262429762, 7)
3.0 (0.99997783214688762, 9)
5.0 (0.99999984846848966, 11)
10.0 (0.99999994150570803, 17)
which are the x values, the estimate of erf(x) and and the number of nodes used in computing
it.
Problem 13
This problem is implemented in the python code c 6 2 p 13.py. When we run that code
we get the following
(0.3564274837735889, 23)
which is the integrals value and the number of nodes used in computing it.
Problem 14
This problem is implemented in the python code c 6 2 p 14.py. When we run that code
we get
0.5 (0.35731366811982967, 4)
1.0 (0.42015837642186937, 5)
2.0 (0.60633780361874334, 5)
which correspond to the arguments to the C(·) function, the value of the integral, and the
number of nodes used in Gauss-Legendre integration in computing the integral.
Problem 15
From the given suggestion we would write
Z π/2
Z 0.01
Z
ln(sin(x))dx =
ln(sin(x))dx +
0
0
63
0.2
ln(sin(x)dx +
0.01
Z
π/2
ln(sin(x)dx
0.2
This first integral we approximate using
Z 0.01
Z 0.01
ln(sin(x))dx ≈
ln(x)dx = x ln(x) − x|00.01 = −0.0560517 .
0
0
This problem is implemented in the python code c 6 2 p 15.py. When we run that code
we get
Integral Parts=
Total Integral=
-0.0560517018599 (-0.46628056164049114, 13) (-0.56646062443647549, 11)
-1.08879288794
Note I’m not really sure we didn’t just split the integral up into two parts (rather than
three). The first from (0, 0.01) and the second from (0.01, π2 ).
Problem 16
See the python code c 6 2 p 16.py where we follow the hint and perform the given integration. When we run that code we get the following output
Pressure Center height=
60.573030
Problem Set 6.3
Problem 1
See the python code c 6 3 p 1.py. When we run that code we get 1.77777777778. Lets
check this result by evaluating this integral exactly. We have
Z +1 Z +1
I=
(1 − x2 )(1 − y 2)dxdy
−1
=
−1
+1 #2 2 2 2
1
1
2
4
x3 = 1 − − −1 +
= 2−
=
= 1.777777 .
x− 3 −1
3
3
3
3
"
The same result.
64
Problem 2
See the python code c 6 3 p 2.py. When we run that code we get 24.0. Lets check this
result by evaluating this integral exactly. We have
3 3
Z 2 Z 3
Z 2
x 2 2
2
I=
x y dxdy =
y
dy
3 0
y=0 x=0
y=0
2
Z 2
27(8)
27 y 3 27
2
=
dy =
= 24 .
=
y
3
3
3
9
y=0
0
The same result.
Problem 3
See the python code c 6 3 p 3.py. When we run that code we get
m=
6, approx=
2.230983, exact=
2.230985, error=
0.000002
Problem 4
See the python code c 6 3 p 4.py. When we run that code we get
m=
3, approx=
1.623391, exact=
1.621139, error=
-0.002252
Problem 5
Our mappings for the four points located at A = (0, 0), B = (2, 0), C = (4, 4) and D = (0, 4)
is given by
x(ξ, η) =
4
X
Nk (ξ, η)xk
k=1
= 0 + 2N2 (ξ, η) + 4N3 (ξ, η) + 0
1
1
=0+2
(1 + ξ)(1 − η) + 4
(1 + ξ)(1 + η)
4
4
1
= (1 + ξ)(3 + η) .
2
65
y(ξ, η) =
4
X
Nk (ξ, η)yk
k=1
= 0 + 0 + 4N3 (ξ, η) + 4N4 (ξ, η)
1
1
(1 + ξ)(1 + η) + 4
(1 − ξ)(1 + η)
=4
4
4
= 2(1 + η) .
With these we have
J(ξ, η) =
"
∂x
∂ξ
∂x
∂η
∂y
∂ξ
∂y
∂η
#
=
1
(3 +
2
1
(1 +
2
η) 0
ξ) 2
.
So that |J(ξ, η)| = 3 + η. Integrating in these transformed coordinates we can evaluate our
integral as
ZZ
I=
xydxdy
Z +1 Z +1
1
=
(1 + ξ)(3 + η)2(1 + η)(3 + η)dξdη
−1
−1 2
+1 Z +1
Z +1
ξ 2 2
(3 + η) (1 + η)dη = 2
(9 + 6η + η 2 )(1 + η)dη
= ξ+ ×
2 −1
−1
−1
Z +1
=2
(9 + 6η + η 2 + 9η + 6η 2 + η 3 )dη dropping odd monomials we get
−1
1
Z +1
7η 3 7
2
=2
(9 + 7η )dη = 4 9η +
=4 9+
= 45.33333 .
3 0
3
−1
See the python code c 6 3 p 5.py where we evaluate this integral numerically and find the
value 45.3333333333 the same as above.
Problem 6
Let the four corners of the given figure be located at A = (0, 0), B = (2, 0), C = (5, 4) and
D = (1, 4). Then our mapping functions for this problem become
x(ξ, η) =
4
X
Nk (ξ, η)xk
k=1
= 0 + 2N2 (ξ, η) + 5N3 (ξ, η) + 1N4 (ξ, η)
1
5
1
= (1 + ξ)(1 − η) + (1 + ξ)(1 + η) + (1 − ξ)(1 + η)
2
4
4
1
= (8 + 4η + 6ξ + 2ξη) .
4
66
y(ξ, η) =
4
X
Nk (ξ, η)yk
k=1
= 0 + 0 + 4N3 (ξ, η) + 4N4 (ξ, η)
1
1
(1 + ξ)(1 + η) + 4
(1 − ξ)(1 + η)
=4
4
4
= 2(1 + η) .
With these we have
J(ξ, η) =
"
∂x
∂ξ
∂x
∂η
∂y
∂ξ
∂y
∂η
#
=
1
(3 +
2
1
(2 +
2
η) 0
ξ) 2
.
So that |J(ξ, η)| = 3 + η. Integrating in these transformed coordinates we can evaluate our
integral as
ZZ
I=
xdxdy
Z +1 Z +1
1
=
(8 + 4η + 6ξ + 2ξη) (3 + η)dξdη
−1
−1 4
Z
Z
1 +1 +1
=
(24 + 12η + 18ξ + 6ξη + 8η + 4η 2 + 6ξη + 2ξη 2 )dξdη
4 −1 −1
Z
Z
1 +1 +1
(24 + 20η + 18ξ + 12ξη + 4η 2 + 2ξη 2 )dξdη
=
4 −1 −1
Z
Z
1 +1 +1
=
(24 + 4η 2 )dξdη
4 −1 −1
Z
Z 1
4
76
2 +1
2
(24 + 4η )dη =
(24 + 4η 2)dη = 24 + =
= 25.33333 .
=
4 −1
3
3
0
See the python code c 6 3 p 6.py where we evaluate this integral numerically and find the
value 25.3333333333 the same as above.
Problem 7
See the python code c 6 3 p 7.py. When we run that code we get 13.5.
Problem 8
See the python code c 6 3 p 8.py. When we run that code we get 24.3.
Problem 9
See the python code c 6 3 p 9.py. When we run that code we get 18.0.
67
Problem 10
See the python code c 6 3 p 10.py. When we run that code we get 7.2.
Problem 11
See the python code c 6 3 p 11.py. When we run that code we get
m=
6, approx=
41.853968
Problem 12
See the python code c 6 3 p 12.py. When we run that code we get
m=
m=
m=
m=
m=
6,
7,
8,
9,
10,
approx=
approx=
approx=
approx=
approx=
0.378838
0.379411
0.379595
0.379547
0.379554
Problem 13
See the python code c 6 3 p 13.py. When we run that code we get -0.0416666666667.
Problem 14
See the python code c 6 3 p 14.py. When we run that code we get
approx=
0.310239, exact=
0.318310, error=
0.008070
Problem 15
See the python code c 6 3 p 15.py. When we run that code we get
m=
2, approx=
-0.194253
68
m=
m=
m=
3, approx=
4, approx=
5, approx=
-0.202803
-0.202641
-0.202642
Problem 16
As suggested we evaluate the integral over the entire domain by breaking it up into the
integral over four triangles. See the python code c 6 3 p 16.py. When we run that code we
get the value 0.133333333333.
69
Chapter 7 (Initial Value Problems)
Problem Set 7.1
Problem 1
Our differential equation for this problem is
y ′ = −4y + x2
with y(0) = 1 .
The Taylor series of order two we will use is
1
y(x + h) ≈ y(x) + y ′(x)h + y ′′ (x)h2 .
2
We use the differential equation to derive the needed derivatives
y ′(x) = −4y(x) + x2
y ′′(x) = −4y ′ (x) + 2x = −4(−4y(x) + x2 ) + 2x = 16y(x) − 4x2 + 2x .
Starting from x = 0 and with a step-size given by h = 0.05 our Taylor series approximation
to the solution at y(0.05) gives
1
y(0.05) = y(0) + y ′(0)(0.05) + y ′′ (0)(0.05)2
2
1
2
= 1 + (−4(1) + 0 )(0.05) + (16(1) − 4(02 ) + 2(0))(0.05)2 = 0.82 .
2
Next we start from x = 0.05 again with h = 0.05 to get
1
y(0.1) = y(0.05) + y ′(0.05)(0.05) + y ′′(0.05)(0.05)2
2
1
= y(0.05) + (−4y(0.05) + 0.052)(0.05) + (16y(0.05) − 4(0.05)2 + 2(0.05))(0.052)
2
= 0.6726375 .
We want to compare this value to that we obtained in Example 7.1 where the exact solution
is y(0.01) = 0.670623.
Problem 2
Our differential equation for this problem is y ′ (x) = −4y(x) + x2 .
Part (a): For the Runge-Kutta method of order two with h = 0.1 when we follow the book’s
70
equations 7.9 we get
k0 = 0.1(−4y(0) + 02 ) = 0.1(−4(1)) = −0.4
" 2 #
1
0.1
k1 = 0.1 −4 y(0) + (−0.4) + 0 +
2
2
= 0.1 −4(1 − 0.2) + 0.052 = −0.39175
y(0.1) = y(0) + k1 = 1 − 0.39175 = 0.60825 .
Part (b): For the Runge-Kutta method of order four with h = 0.1 when we follow the
book’s equations 7.10 we have
k0 = 0.1(−4y(0) + 02 ) = −0.4
" 2 #
0.1
1
= −0.39175
k1 = 0.1 −4 y(0) + (−0.4) + 0 +
2
2
" 2 #
0.1
1
= −0.3358
k2 = 0.1 −4 y(0) + (−0.39175) + 0 +
2
2
k3 = 0.1 −4 (y(0) + k2 ) + (0 + 0.1)2 = −0.26468
1
y(0.1) = y(0) + (k0 + 2k1 + 2k2 + k3 ) = 0.6707033 .
6
Problem 3
For this problem we will use the routine taylor which does fourth order integration (rather
than the requested second order integration) using a Taylor series approximation. To do
the problem as stated we would have to code up a version similar to this routine or do the
calculations by hand. Since we are asked to go from x = 0 to x = 0.5 in steps of h = 0.1 (i.e.
five steps) this seems a bit tedious to do by hand. To code up a version similar to the books
taylor function we would enter the same code as taylor but with the range(4) replaced
with range(2). In either case it seems like just using taylor is the easiest thing to do.
To solve this differential equation with Taylor’s series the needed derivatives are given by
y ′ = sin(y)
y ′′ = cos(y)y ′ = cos(y) sin(y)
y ′′′ = − sin2 (y)y ′ + cos2 (y)y ′ = − sin3 (y) + cos2 (y) sin(y)
y (4) = −3 sin2 (y) cos(y)y ′ − 2 cos(y) sin2 (y)y ′ + cos3 (y)y ′
= −3 sin3 (y) cos(y) − 2 cos(y) sin3 (y) + cos3 (y) sin(y) = −5 sin3 (y) cos(y) + cos3 (y) sin(y) .
This problem is then implemented in the python code c 7 1 p 3.py. When we run that
code we get
x
y[ 0 ]
71
0.0000e+00
5.0000e-01
1.0000e+00
1.4664e+00
Which agrees with the result from the books Example 7.4 where we find that y(0.5) = 1.4664.
Problem 4
Both of the given functions satisfy y(0) = 0 as they must. The function y(x) = 0 satisfies
the differential equation trivially. For the other function we have
1/2
3/2
3 1/2
2
2
′
x =
x1/2 =
y (x) =
3
2
3
3/2 !1/3
2
x
= y 1/3 ,
3
as we were to show.
We can use the routine run kut4 to integrate this differential equation with the two given
initial conditions. When we specify y(0) = 0 the solution y = 0 is numerically computed as
can be seen from the first output from c 7 1 p 4.py
Part (a): y(0) = 0
x
y[ 0 ]
0.0000e+00
0.0000e+00
5.0000e+00
0.0000e+00
If we specify the initial condition y(0) = 10−16 then we get
Part (b): y(0) = 10^(-16)
x
y[ 0 ]
0.0000e+00
1.0000e-16
5.0000e+00
6.0850e+00
we will note that while in this case we see that the numerical value y(5) is equal to (2x/3)3/2
evaluated at x = 5. Thus in this case we are computing the second solution.
Problem 5
Part (a): Here we can write our differential equation as y ′ = exp(sin(x) − y) so F (x, y) =
exp(sin(x) − y).
72
Part (b): For this problem we let
y0 = y
y1 = y ′ .
Then we have
y0′ = y1
y1′ = y ′′ =
2y 2 + xy ′
x
x
= 2y + y ′ = 2y0 + y1 .
y
y
y0
From this our function F (x, y) is given by
F (x, y) =
y1
2y0 + yx0 y1
.
Part (c): For this part our differential equation is given by y (4) = 4y ′′
y0
y1
y2
y3
=y
= y′
= y ′′
= y ′′′ .
p
1 − y 2 so we let
Then we have
y0′ = y1
y1′ = y2
y2′ = y3
y3′ = y (4) = 4y2
From these our function F (x, y) is
p
1 − y0 2 .


y1


y2
.
F (x, y) = 


py3
4y2 1 − y0 2
p
Part (d): For this part our differential equation is given by y ′′ = ± |32y ′x − y 2 |. Let
y0 = y
y1 = y ′ .
Then our differential equations are
y0′ = y1
p
y1′ = y ′′ = ± |32y1 x − y0 2 | .
73
From these our function F (x, y) is given by
y
1
p
.
F (x, y) =
± |32y1x − y0 2 |
We would need to specify which of the signs above to use. To do this we would need to know
what the sign of y ′′ was (hopefully always of one sign) and then select the corresponding sign
in the above expression for F (x, y).
Problem 6
Part (a): We let the vector v have components
v1
v2
v3
v4
=y
= y′
=x
= x′ .
Then get
v1′
v2′
v3′
v4′
= v2
= y ′′ = v3 − 2v1
= v4
= x′′ = v1 − v3 .
Then our system of differential equations is

which is a function of the vector v.

v2
 v3 − 2v1 
d
,
v=


v4
dt
v1 − v3
Part (b): Again we let the vector v have components
v1
v2
v3
v4
=y
= y′
=x
= x′ .
Then get
v1′ = v2
v2′ = y ′′ = −v1 (v22 + v42 )1/4
v3′ = v4
v4′ = x′′ = −v3 (v22 + v42 )1/4 − 32 .
74
Notice that this right-hand-side of these derivatives is a function of the vector v.
Part (c): Not sure how to handle the y ′′2 term in the first equation. I can think of two
ways. The first would be to take a square root to get the equation
p
y ′′ = ± 4x′ − t sin(y) .
To do this we need to know the sign of y ′′ to take in the above expression i.e. to know if y ′′
is positive or negative and hope that it does not change its sign (or that we can determine
when it does change).
Another method might be to take a time derivative of the this equation to get
2y ′′y (3) + sin(y) + t cos(y)y ′ = 4x′′ .
Then we can solve for y (3) to get
y (3) =
1
(4x′′ − sin(y) − t cos(y)y ′) .
2y ′′
In this formulation we would then let the vector v be given by
v1
v2
v3
v4
v5
v6
=y
= y′
= y ′′
=x
= x′
= x′′ .
and our vector differential equation is then given by
 
 

v2
v1
 
 v2  
v3
  1
 

′′
′
 
 
d 
 v3  =  2y′′ (4x − sin(y) − t cos(y)y )  = 



 
v5
dt  v4  
 
 v5  
 
v6
4y ′ −t cos(y)
v6
x
v2
v3
1
(4v6 − sin(v1 ) − t cos(v1 )v2 )
2v3
v5
v6
4v2 −t cos(v1 )
v4
The right-hand-side of the above is a vector function of the form F (t, y).




,



Problem 7
If we let y0 = θ(t) and y1 = θ′ (t) then the system we want to integrate is
d y0
y1
.
=
− sin(y0 )
dt y1
This problem is worked in the python code c 7 1 p 7.py. When we run that code we get
the plot given in Figure 20 (left). Based on visually looking at the first peak the period of
this motion appears to be 6.682 > 2π = 6.283185 (when L = g = 1).
75
500
1.0
400
0.5
300
y(t)
theta(t)
200
0.0
100
0
−0.5
−100
−1.0
0
4
2
time (t)
6
8
−200
0
10
4
2
6
time (t)
8
10
12
14
Figure 20: Left: A plot of the solution θ(t) for Problem 7.1.7. Right: A plot of the solution
y(t) for Problem 7.1.8.
0.35
0.6
0.30
0.5
0.25
0.4
y(t)
y(t)
0.20
0.15
0.3
0.10
0.2
0.05
0.00
0
1
2
3
time (t)
4
5
0.1
0.0
6
0.5
1.0
time (t)
1.5
2.0
Figure 21: Left: A plot of the solution y(t) for Problem 7.1.9. Right: A plot of the solution
y(t) for Problem 7.1.10.
Problem 8
If we let v0 = y(t) and v1 = y ′ (t) then the system we want to integrate is
d v0
v1
.
=
− g − CmD v12
dt v1
This problem is worked in the python code c 7 1 p 8.py. When we run that code we get
the plot given in Figure 20 (right). The location in time where y(t) = 0 is given by linear
interpolation to be
Time when y=0 is=
12.3055352965
76
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
r(t)
y(t)
0.2
0.0
0.0
−0.2
−0.2
−0.4
−0.4
−0.6
−0.6
−0.8
0
2
4
time (t)
6
8
−0.8
0
10
5
10
15
20
time (t)
25
30
35
40
Figure 22: Left: A plot of the solution θ(t) for Problem 7.1.11. Right: A plot of the
solution r(t) for Problem 7.1.12.
Problem 9
If we let v0 = y(t) and v1 = y ′ (t) then the system we want to integrate is
v1
d v0
= P (t)
.
k
−m
v0
dt v1
m
This problem is worked in the python code c 7 1 p 9.py. When we run that code we get the
plot given in Figure 21 (left). Visually from this plot it looks like the maximum displacement
is for y ≈ 2.145.
Problem 10
If we let v0 = y(t) and v1 = y ′ (t) then the system we want to integrate is
d v0
v1
=
.
g(1 − av03 )
dt v1
This problem is worked in the python code c 7 1 p 10.py. When we run that code we get
the plot given in Figure 21 (right). Looking at the plot our period seems to be given by
T ≈ 0.806 with an amplitude of A ≈ 0.49.
Problem 11
If we let v0 = θ(t) and v1 = θ′ (t) then the system we want to integrate is
d v0
v1
=
.
2
− Lg sin(v0 ) + ωL Y cos(v0 ) sin(ωt)
dt v1
This problem is worked in the python code c 7 1 p 11.py. When we run that code we get
the plot given in Figure 22 (left). From the given plot it looks like the largest value off θ is
when t ≈ 5.6 where θ(t) ≈ 0.8.
77
the projectiles path
15
10
y(t)
5
0
−5
−10
0
10
20
30
x(t)
40
50
60
70
Figure 23: A plot of the solution (x(t), y(t)) for Problem 7.1.13. The solution is only valid
when y(t) > 0 since after that point the projectile is on the ground.
Problem 12
This problem is worked in the python code c 7 1 p 12.py. When we run that code we get
the plot given in Figure 22 (right). From that plot I don’t see any points in time where
r(t) > 2 and thus the slider seems to never reach the tip of the rod.
Problem 13
This problem is worked in the python code c 7 1 p 13.py. When we run that code we get
the plot given in Figure 23. From that figure (and using linear interpolation of the computed
solution to the differential equation) we compute
Time when y=0 is=
Projectile range=
3.47323809484
61.8000138411
Problem 14
This problem is worked in the python code c 7 1 p 14.py. When we run that code we get
x
0.0000e+00
5.0000e-01
y[ 0 ]
6.2832e+00
8.3768e+00
y[ 1 ]
0.0000e+00
6.7175e+00
78
1.2
1.0
1.0
0.8
0.6
theta(t)
theta(t)
0.5
0.4
0.2
0.0
0.0
−0.2
−0.5
0.0
0.1
0.2
0.3
t (time)
0.4
0.5
−0.4
0.0
0.6
0.1
0.2
t (time)
0.3
0.4
0.5
Figure 24: Left: A plot of the solution θ(t) for Problem 7.1.15. Right: A plot of the
solution θ(t) for Problem 7.1.16.
for the initial and final values for y0 = θ(t) and y1 = θ′ (t).
Problem 15
This problem is worked in the python code c 7 1 p 15.py. When we run that code we get
the plot given in Figure 24 (left). Running that code we find that the value of r when θ = 0
(for the first time) is given by
R when theta=0 is=
0.616495
Problem 16
This problem is worked in the python code c 7 1 p 16.py. When we run that code we get
the plot given in Figure 24 (right). Running that code we find that the value of r when
θ = 0 (for the first time) is given by
R when theta=0 is=
0.673920
Problem 17
This problem is worked in the python code c 7 1 p 17.py. When we run that code we get
the plot given in Figure 25.
79
0.15
0.10
y(t)
0.05
0.00
−0.05
−0.10
0.0
0.1
0.2
t (time)
0.3
0.4
0.5
Figure 25: A plot of the solution y(t) for Problem 7.1.17. The value of y0 − 4µmg/k is drawn
as a horizontal black line.
Problem 18
This problem is worked in the python code c 7 1 p 18.py. When we run that code we get
the plots given in Figure 26.
Problem 19
3
3
2
2
1
1
y(t)
y(t)
This problem is worked in the python code c 7 1 p 19.py. When we run that code we get
the following
0
0
−1
−1
−2
−2
−3
0
5
10
t (time)
15
−3
0
20
5
10
t (time)
15
20
Figure 26: Left: A plot of the solution y(t) for Problem 7.1.18 Part (a). Right: A plot of
the solution y(t) for Problem 7.1.18 Part (b).
80
x
1.0000e-12
5.0000e+00
y[ 0 ]
1.0000e+00
-1.7760e-01
y[ 1 ]
0.0000e+00
3.2758e-01
which gives the value of J0 (5) and agrees with the book.
Problem 20
Part (a): If we put y = erx into our differential equation (and then divide by it) we get
Thus our solution is given by
r 2 = 16.81 so r = ±4.1 .
y(x) = Ae−4.1x + Be4.1x .
The initial conditions y(0) = 1 and y ′(0) = −4.1 imply that
y(0) = A + B = 1
y ′(0) = −4.1A + 4.1B = −4.1 ,
which have solutions A = 1 and B = 0. Thus the solution that satisfies the initial conditions
is given by
y(x) = e−4.1x .
Part (b): We would anticipate numerical difficulties when solving this problem for if we
have rounding errors in our solution we might end up with the numerical method computing
the e+4.1x solution which as x → ∞ tends to infinity and will swamp the solution e−4.1x
which tends to zero as x → ∞.
Part (c): We would be worried (when using a fixed step size algorithm) if we choose a step
size h larger than the threshold
2
hmax =
.
λmax
In this problem λmax = 4.1 so evaluating the above gives hmax = 0.487805.
In the python code c 7 1 p 20.py we took h = 0.5. When we run that code we get the
plot given in Figure 27. In that plot we see that the approximate solution is diverging from
the true solution. Larger values for h make this problem worse. If we take h small enough
(even with a fixed step size algorithm) we can compute the correct solution. Note that using
an adaptive algorithm (i.e. run kut5) removes these problems. This adds evidence to the
comments made in the book that adaptive algorithms are almost always “better” than non
adaptive counterparts.
Problem 21
This problem is worked in the python code c 7 1 p 21.py. When we run that code we get
the plot given in Figure 28 (left).
81
1.0
approx. sol.
truth
0.8
y(t)
0.6
0.4
0.2
0.0
0
1
2
3
4
t (time)
5
7
6
8
Figure 27: A plot of the approximate solution y(t) (in blue) and the true solution (in green)
for Problem 7.1.20. Note that the approximate solution is diverging from the true solution.
60
i_1(t)
i_2(t)
i_1(t)
i_2(t)
12
40
10
8
20
6
0
4
2
−20
0
−40
0.00
−2
0.01
0.02
t (time)
0.03
0.04
0.05
0.00
0.01
0.02
t (time)
0.03
0.04
0.05
Figure 28: A plot of the approximate solutions i1 (t) and i2 (t) for Left: Problem 7.1.21 and
Right: Problem 7.1.22.
82
Problem 22
This problem is worked in the python code c 7 1 p 22.py. When we run that code we get
the plot given in Figure 28 (right).
Problem Set 7.2
Problem 1
Let y(x) = erx . When we put that into our equation (and divide by it) we get
r 2 + r − 380 = 0 .
This has roots r = 19 and r = −20 so the solution is given by
y(x) = Ae−20x + Be19x .
The initial conditions y(0) = 1 and y ′(0) = −20 become
y(0) = A + B = 1
y ′ (0) = −20A + 19B = −20 .
From which we find A = 1 and B = 0 so the solution that satisfies the initial conditions is
given by
y(x) = e−20x .
We would expect difficulties in solving this problem with a fixed step size algorithm if h was
not take small enough. An adaptive step size algorithm should remove these worries.
Problem 2
Part (a): For the given expression for y(x) we find
y(0) = −0.01 + 10.01 = 10 ,
as it should. Next we compute
y ′ (x) = 0.1 + 10.01(−10)e−10x
= 0.1 − 10(y − 0.x + 0.01) = 0.1 − 10y + x − 0.1 = x − 10y ,
which shows that the y(x) given is a solution to the given initial value problem.
Part (b): Here λmax = 10 so from the discussion in the book we have that for stability we
must take
2
2
=
= 0.2 ,
h<
λmax
10
In practice we would use a smaller fraction of this (say 0.9 times this value).
83
Problem 3
This problem is worked in the python code c 7 2 p 3.py. When we run that code we get
final values for the integrations with each of the different h values as follows
h
0.1
0.25
0.5
x
5.0000e+00
5.0000e+00
5.0000e+00
y[ 0 ]
4.9000e-01
4.9173e-01
2.3457e+12
At x = 5 the true value for y(x) is 0.49. Notice that for h small (less than the 0.2 computed
above) our integration is stable and we get a result that is at least close to the true value of
0.49. Once h gets too large the results we get are garbage due to instability.
Problem 4
This problem is worked in the python code c 7 2 p 4.py. When we run that code we get
final values for the integrations with each of the different h values as follows
h
0.1
0.25
0.5
x
1.0000e+01
1.0000e+01
1.0000e+01
y[ 0 ]
9.9000e-01
9.9000e-01
9.9000e-01
At x = 10 the true value for y(x) is 0.990. Notice that for each value of h the adaptive
algorithm produces the correct result.
Problem 5 & 6
Part (a): If y(t) = ert then the characteristic equation for this differential equation is
r2 +
c
k
r+
= 0.
m
m
Thus given the two roots r1 and r2 to the above quadratic the solution to our differential
equation is
y(x) = Aer1 t + Ber2 t ,
where A and B are chosen to satisfy y(0) = 0.01 and y ′(0) = 0. For the given values of c,
m, and k the roots of our characteristic equation are
r1 = −229.01754251 and r2 = −0.98245749 .
84
Thus we have λmax ≈ 230 so that our step size constraint is that
h < hmax =
2
= 0.008695652 ,
230
if we were to use a fixed step size algorithm on this problem.
Part (b): This problem is worked in the python code c 7 2 p 5.py. When we run that
code we get
x
0.0000e+00
2.0000e-01
y[ 0 ]
1.0000e-02
8.2515e-03
y[ 1 ]
0.0000e+00
-8.1067e-03
Running the code with run kut5 rather than run kut4 (in the same python code) did not
make perceptible differences in the outputs.
Problem 7
We derived analytical solutions to this differential equation in Problem 20 in the previous
section. Using results from that problem we have that the true solutions to each part of this
problem are given by
ya (t) = e−4.1t
yb (t) = 1.00121951e−4.1t − 0.00121951e+4.1t .
Thus as t gets larger the second term in the solution for Part (b) begins to increase in size
and sends the solution to −∞. This problem is worked in the python code c 7 2 p 7.py.
When we run that code it produces the plot given in Figure 29.
Problem 8
This problem is worked in the python code c 7 2 p 8.py. When we run that code we see
a large increase in the value of y as x approaches the value 3.5. Seeing this behavior when
using an adaptive method indicates that this is a real effect and not an artifact of an unstable
method. In fact, this equation is known as a Riccati differential equation and can be solved
analytically. Typically Riccati equations can have singularities inside their domains.
Problem 9
This problem is worked in the python code c 7 2 p 9.py. When we run that code we plot
the numerical solution as requested.
85
1
0
y(t)
−1
−2
−3
−4
−5
0.0
Part (a): approx. sol.
Part (a): truth
Part (b): approx. sol.
Part (b): truth
0.5
1.0
t (time)
1.5
2.0
Figure 29: The true and approximate solutions to the two problems in 7.2.7. Note that the
approximate solution is so close to the true one that the two curves lie on top of each other.
Problem 10
This problem is worked in the python code c 7 2 p 10.py. When we run that code we
plot the numerical solution and the true solution on the same graph. The two curves are
indistinguishable at the resolution displayed.
Problem 11
This problem is worked in the python code c 7 2 p 11.py. When we run that code we plot
the numerical solution as requested.
Problem 12
This problem is worked in the python code c 7 2 p 12.py. When we run that code we plot
the numerical solution as requested. Qualitatively this solution has much different behavior
than the one from the previous problem even though the only thing that is different is the
initial conditions.
86
3 1e7
2
x(t)
1
0
−1
−2
−3
0.0
0.2
0.4
t
0.6
0.8
1.0
Figure 30: The approximate solution to problem in 7.2.16.
Problem 13
This problem is worked in the python code c 7 2 p 13.py. When we run that code we plot
the numerical solution as requested.
Problem 14
This problem is worked in the python code c 7 2 p 14.py. When we run that code we plot
the numerical solution as requested.
Problem 15
This problem is worked in the python code c 7 2 p 15.py. When we run that code we plot
the numerical solution as requested.
87
20
phi'(t)
15
10
5
0.0
0.2
0.4
0.6
t
0.8
1.0
1.2
1.4
Figure 31: The approximate solution to problem in 7.2.17.
Problem 16
This problem is worked in the python code c 7 2 p 16.py. When that code is run it produces
the plot given in Figure 30. From that plot we find (visually) that the amplitude is about
2.9 107 and the period is about 0.59.
Problem 17
This problem is worked in the python code c 7 2 p 17.py. When that code is run it produces
the plot given in Figure 31.
Problem 18
This problem is worked in the python code c 7 2 p 18.py. When that code is run it produces
the plot given in Figure 32.
88
i(t)
2
1
0
−1
−2
0
2
4
t (time)
6
8
10
Figure 32: The approximate solution to problem in 7.2.18.
i_1(t)
i_2(t)
100
50
0
−50
0.00
0.01
0.02
t (time)
0.03
0.04
0.05
Figure 33: The approximate solution to problem in 7.2.19.
89
1.6
i_1(t)
i_2(t)
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.00
0.01
0.02
t (time)
0.03
0.04
0.05
Figure 34: The approximate solution to problem in 7.2.20.
Problem 19
This problem is worked in the python code c 7 2 p 19.py. When that code is run it produces
the plot given in Figure 33.
Problem 20
This problem is worked in the python code c 7 2 p 20.py. When that code is run it produces
the plot given in Figure 34.
Problem 21
This problem is worked in the python code c 7 2 p 21.py. When that code is run it produces
the plot given in Figure 35.
90
y_0(t)
y_1(t)
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0
10
20
t (time)
30
40
50
Figure 35: The approximate solution to problem in 7.2.21.
u(t)
4
2
0
−2
−4
0
c=8.200
c=8.300
2
4
t (time)
6
8
10
Figure 36: The approximate solution to problem in 7.2.22.
91
Problem 22
This problem is worked in the python code c 7 2 p 22.py. When that code is run it produces
the plot given in Figure 36. Notice that the final value where u(t) ends is quite different
between the two runs. Numerical solutions for the two different values of c are given by
c=
8.2
x
0.0000e+00
1.0000e+01
c= 8.3
x
0.0000e+00
1.0000e+01
y[ 0 ]
0.0000e+00
1.5271e+00
y[ 1 ]
1.0000e+00
6.6165e-01
y[ 2 ]
2.0000e+00
8.5126e+00
y[ 0 ]
0.0000e+00
-4.7495e+00
y[ 1 ]
1.0000e+00
-4.9843e+00
y[ 2 ]
2.0000e+00
8.9369e+00
92
Chapter 8 (Two-Point Boundary Value Problems)
Problem Set 8.1
Problem 8
If the solution goes from (0, 0) to (π/2, 1) in a straight line it would have a slope
1−0
2
= = 0.6366198 .
π
π
−0
2
Based on this we might take as initial slopes that (hopefully) bracket the true value of y ′ (0) at
the values 0.5 and 1.0. See the python code c 8 1 p 8.py where we implement this problem.
Running that code produces a plot of the solution function y(x).
Problem 9
If the solution goes from (0.01, 0) to (2, −11) in a straight line it would have a slope
−1 − 0
= −0.5025126 .
2 − 0.01
Based on this we might take as initial slopes that (hopefully) bracket the true value of y ′(0)
at the values -0.2 and -1.0. See the python code c 8 1 p 9.py where we implement this
problem. Running that code produces a plot of the solution function y(x).
Problem 10
From the differential equation we have
y ′′ = − sin(y) − 1 .
If 0 < y < π over the domain 0 < x < π then sin(y) > 0 and y ′′ < −1 < 0. A negative
second derivative means that our function y(x) is concave downward function that starts at
the point (0, 0) and ends at the point (π, 0). Thus we might hypothesize that two slopes that
1
(a small value) and +10 (a large value). See the
bracket the true value of y ′(0) would be 10
python code c 8 1 p 10.py where we implement this problem. Running that code produces
a plot of the solution function y(x).
Problem 11
For the function to end with y ′(2) = 0 means that when x = 2 our slope of y(x) is flat. Our
function y(x) starts at the point (0.01, 1) but we don’t know the level (or y value) at x = 2
93
since it could be above or below the value y = 1 which is the value on the left-hand-side of
the domain. Thus the initial slope at x = 0.01 could be positive or negative depending. Thus
we might pick two the slopes that hopefully bracket the true value of y ′(0.01) to be −10 and
+10. See the python code c 8 1 p 11.py where we implement this problem. Running that
code produces a plot of the solution function y(x) and the first derivative y ′(x). Using these
two initial slopes gave an initial value of y ′(0.01) ≈ −22.58 which was outside our [−10, +10]
initial bracket on the slopes. Thus we change our bracket on the initial slope to the range
[−100, 0]. As a final note since the solution changes rapidly for x near zero it seems safer
to use an adaptive numerical method like run kut5 rather than the fixed step size method
run kut4.
Problem 12
For the function to start at (0, 1) and then end at the point (∞, 0) means that the function
must decrease as x increases and thus we expect y ′(0) < 0. Thus we can try to bracket our
initial slope by assuming it will be between [−10, 0]. See the python code c 8 1 p 12.py
where we implement this problem. In that code I initially took β = 10 and was able to get
a valid solution (the initial slope y ′(0) was within the initial bracket). The output from that
run gave the values (partial)
beta=
10.000000
x
0.0000e+00
3.8765e-01
1.0390e+00
...
7.3523e+00
9.4134e+00
1.0000e+01
y[ 0 ]
1.0000e+00
7.6180e-01
4.4481e-01
y[ 1 ]
-6.3455e-01
-5.7874e-01
-3.9073e-01
9.0377e-04
7.9901e-05
3.5662e-10
-9.1269e-04
-1.5150e-04
-1.2871e-04
Running the code again with β = 15 gave the following (partial)
beta=
15.000000
x
0.0000e+00
3.8765e-01
1.0390e+00
...
9.4134e+00
1.2364e+01
1.5000e+01
y[ 0 ]
1.0000e+00
7.6180e-01
4.4481e-01
y[ 1 ]
-6.3455e-01
-5.7874e-01
-3.9073e-01
1.1572e-04
6.4306e-06
4.5482e-06
-1.1568e-04
-5.8454e-06
3.5657e-06
Since the selected value for y ′(0) did not change with this increased value of β we can be
satisfied that we have found the correct solution.
94
Problem 13
Based on the fact that the solution y(x) goes from (1, 0) to (2, 1) we assume that y ′ (0) > 0
and assume that its value will be between [0, 10.0]. Note that solving this boundary value
problem using the shooting method is still a one dimensional root finding problem in that we
specify a value for y ′(1) that finds the root of y(2) − 1. See the python code c 8 1 p 13.py
where we implement this problem. Running this code gives that the true value of y ′ (0) is
0.9011 showing that the initial bracket of y ′(0) was perhaps too large.
Problem 14
See the python code c 8 1 p 14.py where we implement this problem. Running this code
plots the solution y(x).
Problem 15
See the python code c 8 1 p 15.py where we implement this problem. Running this code
plots the solution y(x).
Problem 17
For this problem we have to specify values for y ′ (0) and y ′′′(0) thus we need to use a
two dimensional Newton solver to find these two initial conditions. See the python code
c 8 1 p 17.py where we implement this problem. Running this code plots the solution y(x).
Problem 18
For this problem we have to specify values for y ′′ (0) and y ′′′ (0) thus we need to use a
two dimensional Newton solver to find these two initial conditions. See the python code
c 8 1 p 18.py where we implement this problem. Running this code plots the solution y(x).
Problem 19
For this problem we have to specify the value for θ and v0 . See the python code c 8 1 p 19.py
where we implement this problem. Running this code plots the solution (x(t), y(t)) and we
find that
95
u[0] (deg)=
3.362746; u[1]= 854.961865
x
y[ 0 ]
y[ 1 ]
0.0000e+00
0.0000e+00
0.0000e+00
1.0000e+01
8.0000e+03
-7.5512e-06
y[ 2 ]
8.5349e+02
7.5089e+02
y[ 3 ]
5.0150e+01
-4.8051e+01
For the angle on inclination, the launch velocity v0 , and the solution (x(t), y(t)) at the
beginning and end of the trajectory.
Problem 20
See the python code c 8 1 p 20.py where we implement this problem. Running this code
plots the solution y(x).
Problem 21
For this problem we have to specify the value for y ′′(0). To specify the right-hand boundary
we let β = 10 and solve the differential equation on [0, β]. We then look to increase β and see
if the results change. See the python code c 8 1 p 21.py where we implement this problem.
Running this code plots the solution y(x). In that code we find that y ′′ (0) = 1.3282.
96
Chapter 9 (Symmetric Matrix Eigenvalue Problems)
Notes on Jacobi Diagonalization
Here are some of the algebraic steps in the discussion on Jacobi diagonalization. We start
with the quadratic equation t2 +2φt−1 = 0. Which using the quadratic formula has solutions
p
p
−2φ ± 4φ2 − 4(−1)
t=
= −φ ± φ2 + 1 .
2
The more stable root is the one for which |t| ≤ 1 and so we can enforce that by taking
p
t = sgn(φ)(−|φ| + φ2 + 1)
!
p
2+1
p
|φ|
+
φ
p
= sgn(φ)(−|φ| + φ2 + 1)
|φ| + φ2 + 1
!
−φ2 + (φ2 + 1)
sgn(φ)
p
p
= sgn(φ)
=
,
|φ| + φ2 + 1
|φ| + φ2 + 1
which is the books Equation 9.16a.
Problem Set 9.1
Problem 1
Here B is positive definite so we can use the python routine stdForm or just note that since B
is diagonal we can write B = LLT with L = B 1/2 (with the square root taken elementwise).
Then the standard form coefficient matrix H is given by H = L−1 A(L−1 )T which here equals
[[ 1.75
[ 0.5
[ 0.25
0.5
1.
1.
0.25]
1. ]
2. ]]
This gives us the standard form of the eigenvalue problem Hz = λz. Note that the variables
x and z are related via x = (L−1 )T z. From the specific form of L in this case this is
1
xi = √ zi
Bii
for i = 1, 2, . . . , n − 1, n .
This problem is worked in the python code c 9 1 p 1.py.
97
Problem 2
For this problem we will use the python code stdForm on the given matrices A and B to get
the matrices H and T such that the new problem is in standard form where Hz = λz with
x = T z. This is implemented in the python code c 9 1 p 2.py. When we run that code we
get
H=
[[ 2.
0.57735027
0.81649658]
[ 0.57735027
2.66666667
2.3570226 ]
[ 0.81649658
2.3570226
13.33333333]]
T=
[[ 0.70710678 0.40824829 0.57735027]
[ 0.
0.81649658 1.15470054]
[ 0.
0.
1.73205081]]
Problem 3
We will solve this problem with the inverse power method with eigenvalue shifting. This
means that we first “shift” A to get A∗ where
A∗ = A − 2.5I .


1
We initialize the vector v to a random initial vector say  0 . Then we repeat the following
0
steps
• Solve
A∗ z = v .
(13)
for the vector z.
• Compute the magnitude of z or |z| which is our current best guess at
of the smallest eigenvalue.
• Set v =
λ1 .
z
|z|
1
λ1
the reciprocal
this is our current best guess at the eigenvector associated with eigenvalue
• Repeat this sequence of steps above again starting by solving Equation 13 using this
new value of v.
We could do two iterations of this sequence “by hand” or we could use the python code
inversePower. We follow this later approach in the python code c 9 1 p 3.py. When we
run that code we get
98
2.46192411534 [ 0.59371536
0.77601663 -0.21283858]
We would then need to apply the transformation T z to get the eigenvector to the original
problem which we find to be
[ 0.61374521
0.3878501
-0.36864723]
Problem 4
Given the matrix S we want to compute the eigenvalues of S. Since S is a “small” matrix
(less than 20 × 20) we will use the Jacobi method to compute all eigenvalues of S. We
implement this in the python code c 9 1 p 4.py. When we run that code we get
eigenvalues=
[ 73.15341562
80.
eigenvectors=
[[ 0.61541221 0.
[ 0.78820544 0.
[ 0.
1.
196.84658438]
0.78820544]
-0.61541221]
0.
]]
Problem 5
Let θi (t) = Ai sin(ωt) then our equations become
kL(A2 − A1 ) − mgA1 = −mLA1 ω 2
−kL(A2 − A1 ) − 2mgA2 = −2mLA2 ω 2 .
In matrix form these are
A1
mL
0
A1
−mg − kL
kL
2
.
= −ω
A2
0 2mL
A2
kL
−kL − 2mg
To find the eigenvalues and eigenvectors for this problem we will use the stdForm and jacobi.
The way the above is written we cant use stdForm directly since the matrix on the righthand-side is not positive definite. We can multiply both sides by −1 and get a right-hand-side
that is positive definite. We implement all of this in the python code c 9 1 p 5.py. When
we run that code we get
eigenvalues= [ 13.07553333
eigenvectors=
[[ 0.70710678 0.89442719]
[ 0.70710678 -0.4472136 ]]
133.07553333]
The circular frequencies ω would be the square root of the eigenvalues above.
99
Problem 6
Let ik = Ak sin(ωt) then our three equations become
3A1 − A2 − A3 = +LCω 2 A1
−A1 + A2 = +LCω 2 A2
−A1 + A3 = +LCω 2 A3 .
In matrix form these looks like





3 −1 −1
A1
A1
 −1 1
0   A2  = LCω 2  A2  .
−1 0
1
A3
A3
We implement this in the python code c 9 1 p 6.py. When we run that code we get
eigenvalues= [ 0.26794919 1.
eigenvectors=
[[ 4.59700843e-01
2.16337418e-12
[ 6.27963030e-01
7.07106781e-01
[ 6.27963030e-01 -7.07106781e-01
3.73205081]
8.88073834e-01]
-3.25057584e-01]
-3.25057584e-01]]
Once we have found our eigenvalues λi the circular frequencies ωi are the solutions to
λi = LCωi2 .
Problem 7
We want to take the given matrix A and annihilate the two elements A14 and A41 by using
a Jacobi rotation. To do this we will follow the steps on Page 331 of the book with k = 1
and l = 4. Do this we get
A11 − A44
(4 − 4)
=−
=0
2A14
2(1)
t = −1 following ”we choose the plus sign if φ > 0 and the minus sign if φ ≤ 0”
1
1
=√
c= √
1 + t2
2
1
s = tc = − √
2
√
1
−1/ 2
s
√ = −√
.
=
τ=
1+c
1 + 1/ 2
2+1
φ = cot(2θ) = −
100
Now we modify the elements in the upper half of A using the books equation 9.18. These
equations are
A∗11 = A11 − (−1)A14 = 4 − (−1)1 = 5
A∗44 = A44 + (−1)A14 = 4 + (−1)1 = 3
A∗14 = A∗41 = 0
1
1
∗
∗
A1i
A1i = Ai1 = A1i + √ A4i − √
2
2+1
1
1
∗
∗
A4i = Ai4 = A4i − √ A1i + √
A4i
2
2+1
for i 6= 1, i 6= 4
(14)
for i 6= 1, i 6= 4 .
(15)
Evaluating Equation 14 for i = 2 gives
A∗12
1
1
A12
=
= A12 + √ A42 − √
2
2+1
1
1
1
1
= −1 + √ 0 − √
(−1) = −1 + √ √
.
2
2+1
2 2+1
A∗21
Evaluating Equation 14 for i = 3 gives
A∗13
1
1
A13
=
= A13 + √ A43 − √
2
2+1
√
1
2
1
(0) = + √ = 2 .
=0+ √ 2− √
2
2+1
2
A∗31
Evaluating Equation 15 for i = 2 gives
A∗42
1
1
A42
=
= A42 − √ A12 + √
2
2+1
1
1
1
(0) = √ .
= 0 − √ −1 + √
2
2+1
2
A∗24
Evaluating Equation 15 for i = 3 gives
A∗43
1
1
=
= A43 − √ A13 + √
A43
2
2+1
√
1
2
1
(2) = 2 − √
.
=2− √ 0+ √
2
2+1
2+1
A∗34
With all of its components specified we find that A∗ is given by

√
5
−1 + √2(√12+1)
2
0

√1
6
−2
 −1 + √2(√12+1)
2√
√
A∗ = 
2

√
2
−2
3
2 − 2+1

√
2
√1
2 − √2+1
3
0
2
101






Problem 8
See the python code c 9 1 p 8.py for an implementation of this problem. When we run
that code we get the following output
eigenvalues= [-1.38719961 2.69161093 6.69558869]
eigenvectors=
[[ 0.21138388 0.76356158 -0.61015618]
[-0.51841476 0.61680524 0.59228156]
[ 0.82859097 0.19111519 0.52622428]]
Problem 9
See the python code c 9 1 p 9.py for an implementation of this problem. When we run
that code we get the following output
eigenvalues= [ 1.38196601 2.45861873 3.61803399 8.54138127]
eigenvectors=
[[ 0.37174803 -0.5395366
0.60150096 0.45705607]
[ 0.60150096 -0.45705607 -0.37174803 -0.5395366 ]
[ 0.60150096 0.45705607 -0.37174803 0.5395366 ]
[ 0.37174803 0.5395366
0.60150096 -0.45705607]]
Problem 10
We will implement the power method following the python code in for the inverse power
algorithm. This is done in the python code powerMethod.py which is called from the code
c 9 1 p 10.py. When we run this later code we get (partial result)
largest eigenvalue (powerMethod)= 8.54138126515
corresponding eigenvector (powerMethod)=
[ 0.45705607 -0.5395366
0.53953661 -0.45705607]
Whether you convert to this eigenvector or the negative of it depends on the random initial
guess chosen in powerMethod.py. Notice that numerically these match the numbers found
in the previous problem (as they should).
102
Problem 11
This problem is implemented in the python code c 9 1 p 11.py. There I used the inversePower
method from the book with a default argument s = 0 to find the smallest eigenvalue of the
matrix A. When we run the given code we get
1.38196601125 [ 0.37174751
0.60150051
0.6015014
0.37174855]
for the smallest eigenvalue and its associated eigenvector. Notice that these results match
the numbers found in problem 9 (as they should).
Problem 12
For this problem we can follow the outline of the code given in Part (2) in Example 9.3 from
the book. This problem is worked in python code c 9 1 p 12.py. When we run that code
we get
eigenvalues= [ 2.92765173
9.90287536
eigenvectors=
[[ 0.98102277 -0.18709892 0.32245473]
[-0.18761435 -0.46142271 0.78544019]
[-0.04894062 0.86722724 0.52830546]]
25.59804434]
Problem 13
For this problem we have the generalized eigenvalue problem Ax = λBx and inversePower
code given for the book is for the eigenvalue problem in standard form Ax = λx. Thus we
need to use stdForm to transform our problem into standard form and then use inversePower
on it. We do this in the python code c 9 1 p 13.py. When we run that code we get
2.92765172791 [ 0.98102277 -0.18761435 -0.04894063]
For the eigenvalue and associated eigenvectors. Notice that these results match the numbers
found in problem 12 (as they should).
Problem 14
This problem is worked in python code c 9 1 p 14.py. When we run that code we get the
requested eigenvalues and eigenvectors.
103
Problem 15
For this problem since B is not positive definite we cannot simply use the routine stdForm
because it uses the Choleski decomposition to decompose B into LLT . Instead we can
use the LU decomposition to decompose B and convert our problem to the standard form
eigenvalue problem Ax = λx that we can solve using jacobi (or other technique). To do
this we write B = LU then our eigenvalue problem of Ax = λBx becomes
Ax = λLUx .
Next we let v = Ux then x = U −1 v and the above becomes
AU −1 v = λLv .
Next we multiply on the left by L−1 to get
L−1 AU −1 v = λv ,
which is the standard form for the eigenvalue problem. Note that if we want to use some
of the methods developed in this chapter of the book (like jacobi) we would need to check
that the leading matrix L−1 AU −1 is symmetric. To get the eigenvectors of Ax = λBx note
that if v is an eigenvector of the coefficient matrix L−1 AU −1 . Then U −1 v is an eigenvector
of the original problem.
Operationally, to solve this problem we could use LUdecomp to perform the LU decomposition
of B and then form the matrix L−1 AU −1 and call the jacobi routine on this matrix (assuming
it is symmetric) to get the eigenvalues. The code in the book performs the LU decomposition
“in-place” which means that A is overwritten as the algorithm works. This makes extracting
U and L more difficult and perhaps more tedious to code. Thus in solving this problem I
used the numpy linear algebra codes (rather than the codes developed for this chapter) to
perform the sequence of steps suggested above. This sequence of steps could be performed
using only the codes from this book but the subsequent coding would be more tedious.
This problem is worked in the python code c 9 1 p 15.py. When we run that code we get
the following
eigenvalues=
[-7.29608125
eigenvectors=
[[-0.90129973
[-0.23934404
[ 0.35314424
[ 0.07524866
0.10397837
0.92887446
1.31712064]
-0.32873623 -0.42200594 0.06652392]
-0.633369
-0.44310166 0.82256316]
-0.62620499 -0.28241205 -0.43554467]
-0.31407566 0.73879316 -0.35953477]]
Problem 16
Part (a): From the given problem description we have an eigenvalue problem of the form
Ax = λBx. The form of the matrix A is specified in the problem statement. The matrix B
104
is the identity matrix with the value at the (n, n)th position 12 . Given the form of the matrix
B we know that for L to be such that B = LLT we must have L the identity matrix except
in the (n, n)th position which would have the value √12 . With this decomposition of B when
we form the H matrix using H = L−1 A(L−1 )T for the given L we have that
Hij = Aij ,
for all i and j except in the last row and column of H where we have
√
Hn,j = 2An,j for 1 ≤ j ≤ n − 1 ,
and
Hi,n =
√
2Ai,n
for 1 ≤ i ≤ n − 1 ,
and finally Hn,n = 2An,n . Then P is the matrix such that y = P z = (L−1 )T z or in terms of
the components of z and y we have
yi = zi for 1 ≤ i ≤ n − 1 and
√
yn = 2zn .
Part (b): For this part we will use the Jacobi method on the matrix H to find the eigenvalues
and eigenvectors of H. Denote these eigenvectors by z. Then P z gives the eigenvectors of the
original problem and the eigenvalues of the standard form problem are also the eigenvalues
of the original problem. Once we have the numerical values for λi we set ωi equal to
EI
2
ωi =
n4 λi ,
γL4
which we can take the square root of to get values of ωi . This problem is worked in the
python code c 9 1 p 16.py. When we run that code we get the following
lowest two eigenvalues lambda_i=
corresponding eigenvectors=
[[ 0.01001677 0.0590165 ]
[ 0.03731826 0.17920964]
[ 0.07916782 0.30635888]
[ 0.13287418 0.39424767]
[ 0.19584227 0.41034224]
[ 0.2656386
0.33971697]
[ 0.34006771 0.18577323]
[ 0.41725709 -0.03291499]
[ 0.49574761 -0.28947655]
[ 0.57458736 -0.55851038]]
scaled two values of omega_i=
[ 3.48659656 21.13353825]
[ 0.00121564
0.04466264]
Running the above code with n = 100 gives (partial output)
105
scaled two values of omega_i=
[ 3.51571779 22.02503253]
this first number matches the one given in the book quite well.
Problem 17
This is a problem of the form Au = λBu where A is a tridiagonal matrix and B is a 10 × 10
diagonal matrix with the diagonal specified as in the problem. Once we have found the
eigenvalues λi to that problem the buckling loads Pi are given by
2
L
Pi
λi =
.
EI0 20
To solve this problem we used the jacobi, inversePower, and inversePower3 routines.
This problem is worked in the python code c 9 1 p 17.py. When we run that code we get
the following output (partial)
inversePower lambda= 0.114316054458
inversePower eigenvector= [ 0.14350833 0.27061136 0.36677916 0.42101822 0.42712814
0.40068632 0.35134206 0.28191578 0.19637575 0.09961127]
inversePower3 lambda= 0.114316054458
inversePower3 eigenvector= [ 0.14350828 0.27061128 0.36677909 0.42101818 0.42712814
0.40068636 0.35134212 0.28191585 0.19637581 0.0996113 ]
Problem 18
P
This is the problem to find the eigensystem for Aθ = kL
Bθ, where the matrix B is given by


1 1 1
B= 0 1 1 .
0 0 1
For this matrix B we find that
B −1


1 −1 0
=  0 1 −1  .
0 0
1
Then the eigensystem we seek is equivalent to that
P
θ. Next compute B −1 A where we find
kL


1 −1 0
6 5
−1



3 3
B A = 0 1 −1
0 0
1
1 1
106
of the standard form problem B −1 Aθ =
 

3
3 2 1
2 = 2 2 1 ,
1
1 1 1
which is a symmetric matrix and we can use the methods of this chapter on it. To solve
this problem we will use the inverse power method to find the smallest eigenvalue and the
corresponding eigenvector. Another method would just be to use the jacobi method to
compute all the eigenvalues and eigenvectors and select the smallest one. As the former
approach seems easier that’s what I’ll do. Towards this end this problem is worked in the
python code c 9 1 p 18.py. When we run that code we get the following
(0.30797852836990414, array([-0.32798528,
0.73697623, -0.59100905]))
for the eigenvalue and corresponding eigenvector.
Problem 19
For this problem we let ui = Ai sin(ωt) then our system of differential equations become
m 2
ω A1
k
3m 2
A1 − 2A2 + A3 = −
ω A2
k
2m 2
ω A3 .
A2 − 2A3 = −
k
−2A1 + A2 = −
As a matrix system the above looks like






1 0 0
A1
−2 1
0
A1
 1 −2 1   A2  = − m ω 2  0 3 0   A2  .
k
0 0 2
A3
0
1 −2
A3
Once we find find the eigenvalues λi for this system the circular frequencies ωi are given by
λi = −
m 2
ω .
k i
We solve this problem following the steps in using the jacobi routine. This problem is
worked in the python code c 9 1 p 19.py. When we run that code we get the following
eigenvalues (lambda_i)= [-2.23292464 -1.18092053 -0.2528215 ]
eigenvectors x (to Ax = lambda B x)=
[[ 0.96983375 -0.38362421 0.42955167]
[-0.22589817 -0.31421871 0.75050344]
[ 0.0916107
0.86838878 0.502225 ]]
107
Problem 20
For this problem let ik = Ak sin(ωt) then the differential equations become
2
1
A1 + (A1 − A2 ) = 0
C
C
2
3
−Lω 2 A2 + (A2 − A1 ) + (A2 − A3 ) = 0
C
C
4
3
−Lω 2 A3 + (A3 − A2 ) + (A3 − A4 ) = 0
C
C
5
4
−Lω 2 A4 + (A4 − A3 ) + A4 = 0 .
C
C
−Lω 2 A1 +
In matrix form this is
 
3 −2 0
0
A1



A2   −2 5 −3 0
LCω 2 
 A3  =  0 −3 7 −4
0
0 −4 9
A4


A1
  A2 


  A3  .
A4

To find the eigensystem for this problem it might be easier to use the Jacobi method to
compute all the eigenvalues and eigenvectors directly. Other computational method exist
that take advantage of the tridiagonal structure of the leading coefficient matrix. This
problem is worked in the python code c 9 1 p 20.py. When we run that code we get the
following
eigenvalues=
eigenvectors=
[[ 0.59534188
[ 0.62357255
[ 0.45424838
[ 0.22446322
[
0.90516148
3.38984838
7.06445882
12.64053133]
0.71883111 -0.35460092 -0.05575017]
-0.14011757 0.72063041 0.26873064]
-0.55442425 -0.25950332 -0.64724817]
-0.39530072 -0.53629099 0.71115791]]
Problem 21
For this problem let ik = Ak sin(ωt) then our differential equations become
−LCω 2 (A1 + A1 − A2 ) + A1
−LCω 2 (A2 − A1 + A2 − A3 ) + 2A2
−LCω 2 (A3 − A2 + A3 − A4 ) + 3A3
−LCω 2 (A4 − A3 + A4 ) + 4A4
=0
=0
=0
= 0.
If we simplify some we get
−LCω 2 (2A1 − A2 ) = −A1
−LCω 2 (−A1 + 2A2 − A3 ) = −2A1
−LCω 2 (−A2 + 2A3 − A4 ) = −3A3
−LCω 2 (−A3 + 2A4 ) = −4A3 .
108
As a matrix system this is

2 −1 0
0

−1 2 −1 0
LCω 2 
 0 −1 2 −1
0
0 −1 2
 
1
A1
  A2   0
 

  A3  =  0
0
A4

0
2
0
0
0
0
3
0

A1
0


0   A2
0   A3
A4
4


.

Note that this system is tridiagonal so we could use specialized method that use this structure
to compute the eigensystem of the above, but it might be easier to compute the eigensystem
using the jacobi method which is a more general method. This problem is worked in the
python code c 9 1 p 21.py. When we run that code we get the following
eigenvalues=
eigenvectors=
[[ 0.29523669
[ 0.54736531
[ 0.63965025
[ 0.45174557
[ 0.14601189
0.55385739
1.07879007
2.38800732]
-0.44755364 -0.54195207 0.92984018]
-0.64722638 -0.49925163 -0.36078479]
-0.12995689 0.62062421 0.07170369]
0.60324536 -0.26806965 -0.00949463]]
Problem 22
This problem is implemented in the python code LRAlgorithm.py which is called from the
code c 9 1 p 22.py. When we run that we get the following
LR output= [0.68744042202776923, 2.3856315668593071, 7.9269280111129223]
eigenvalues (jacobi)= [ 0.68744042 2.38563157 7.92692801]
Which compares the results between our LR algorithm implementation and the Jacobi algorithm.
Problem Set 9.2
Problem 1
Recall that from Gershgoren theorem the smallest eigenvalue λmin is bounded below by
λmin ≥ min(ai − ri ) ,
i
(16)
and the largest eigenvalue λmax is bounded above by
λmax ≤ max(ai + ri ) .
i
109
(17)
In these expressions ai is the value of the diagonal element in row i and ri is the sum of the
absolute values of the off diagonal elements in the ith row. Using these we can answer each
part of this problem.
Part (a): For this matrix we have
a1 = 10 r1 = 5 a1 − r1 = 5 a1 + r1 = 15
a2 = 2 r2 = 7 a2 − r2 = −5 a2 + r2 = 9
a3 = 6 r3 = 4 a3 − r3 = −2 a3 + r3 = 10 .
Thus we have
λmin ≥ min(ai − ri ) = −5
i
λmax ≤ max(ai + ri ) = 15 .
i
Part (b): For this matrix we have
a1 = 4 r1 = 4 a1 − r1 = 0 a1 + r1 = 8
a2 = 5 r2 = 5 a2 − r2 = 0 a2 + r2 = 10
a3 = 4 r3 = 5 a3 − r3 = −1 a3 + r3 = 9 .
Thus we have
λmin ≥ min(ai − ri ) = −1
i
λmax ≤ max(ai + ri ) = 10 .
i
Problem 2
For this problem we want to compute the Sturm sequence for this matrix with x = 4 and
x = 2. This is a tridiagonal matrix with the given values on the diagonal and sub/super
diagonal. For each call to sturmSeq we count the number of sign changes in the sequence
P0 (x), P1 (x), · · · , Pn−1 (x), Pn (x) which tells us the number of eigenvalues of this tridiagonal
matrix that are less than x. This is done worked in the python code c 9 2 p 2.py. When
we run that code we get the following
sturmSeq(d,c,2)=
sturmSeq(d,c,4)=
[ 1. 3. 2. 1. -5.]
[ 1.
1. -4. -1. 15.]
From this we see that there is one sign change in the Sturm sequence P0 (2), P1 (2), · · · , Pn−1 (2), Pn (2)
and two sign changes in the sequence P0 (4), P1 (4), · · · , Pn−1(4), Pn (4). Thus we know that
there is one eigenvalue less than x = 2 and two less than x = 4 indicating that there is one
eigenvalue in [2, 4].
110
Problem 3
For this problem we can use the routine lamRange to compute brackets around each of the
eigenvalues. This problem is worked in the python code c 9 2 p 3.py. When we run that
code we get the following brackets on our eigenvalues of A
lamRange(d,c,3)=
[ 2.
3.3125
4.625
5.5
]
Problem 4
For this problem we can use the routine lamRange to compute brackets around each of the
eigenvalues. This problem is worked in the python code c 9 2 p 4.py. When we run that
code we get the following brackets on our eigenvalues of A
lamRange(d,c,3)=
[
5.
6.453125
7.90625
10.8125
]
Problem 5
For this problem we can use the routine lamRange to compute brackets around each of the
eigenvalues. This problem is worked in the python code c 9 2 p 5.py. When we run that
code we get the following brackets on our eigenvalues of A
lamRange(d,c,4)=
[ 0.
0.703125
1.40625
2.8125
3.75
]
Problem 6
For this problem we will call the routine householder and then follow the code on Page 364365 of the book. This problem is worked in the python code c 9 2 p 6.py. When we run
that code we get the following
diagonal d=
[ 12.
14.04
9.96]
subdiagonal c=
[-5.
-3.72]
transformation matrix P=
[[ 1.
0.
0. ]
[ 0. -0.8 -0.6]
[ 0. -0.6 0.8]]
111
Problem 7
For this problem we will call the routine householder and then follow the code on Page 364365 of the book. This problem is worked in the python code c 9 2 p 7.py. When we run
that code we get the following
diagonal d=
[ 4.
6.66666667
subdiagonal c=
[ 2.44948974 -1.74801475
transformation matrix P=
[[ 1.
0.
[ 0.
-0.81649658
[ 0.
0.40824829
[ 0.
-0.40824829
2.93333333
2.4
]
0.73484692]
0.
0.
]
-0.54494926 -0.19069252]
-0.77849894 0.47673129]
0.31139958 0.85811633]]
Problem 8
Note that this is a tridiagonal matrix and as such we can use the routine eigenvals3 to
compute these with N = 5 (the dimension of the matrix A). This problem is worked in the
python code c 9 2 p 8.py. When we run that code we get the following
[
1.43562719
2.93780551
4.15153739
7.44731309
11.02771683]
for the eigenvalues of the matrix A.
Problem 9
Note that this is a tridiagonal matrix so we can use eigenvals3 with N = 2 to find the two
smallest eigenvalues. This problem is worked in the python code c 9 2 p 9.py. When we
run that code we get the following
[ 0.88282382
3.63182668]
Problem 10
We will use householder transformations to transform this matrix to a tridiagonal form.
Once in this form we can use the routine eigenvals3 with N = 3 to determine the three
112
smallest eigenvalues. Once we have the eigenvalues we can use inversePower3 to determine
the associated eigenvector (For simplicity I only do this for the smallest eigenvalue since
the procedure would be repeated for all other eigenvalues). This problem is worked in the
python code c 9 2 p 10.py. When we run that code we get the following (partial)
lowest_three_evalues= [ 3.37684507 4.75931985
eigenvector of lowest eigenvalue= [ 0.69038451
-0.02413926
6.0368161 ]
0.71514996 0.0950872
0.02568968 -0.04056509]
Problem 11
We will use householder transformations to transform this to tridiagonal form. In this form
we can use the routine eigenvals3 with N = 2 to determine the two smallest eigenvalues.
Note we call eigenvals3 use the diagonal d and sub/super diagonal c that are output from
the routine householder. This problem is worked in the python code c 9 2 p 11.py. When
we run that code we get the following
smallest_two_evalues=
[
1.81260083e-06
1.25707571e-05]
Problem 13
2
Since this eigenvalue problem is of the form Ay = mω
By, to start we want to multiply
k
on the left-hand-side by the inverse of


1 0 0
B= 0 3 0 ,
0 0 2
to get our eigenvalue problem in standard form. The coefficient matrix of the standard form
eigenvalue problem or B −1 A is


1 −1 0
 −1 2/3 −1  .
0 −1 1
Note that this is a tridiagonal system. Thus we can use the routine eigenvals3 with N = 3
to get its three eigenvalues. Once we have these we can use inversePower3 to compute
the three eigenvectors by supplying the output form eigenvals3 as initial guesses at the
eigenvalues. This problem is worked in the python code c 9 2 p 13.py. When we run that
code we get the following
evalues= [-0.59066729 1.
eigenvector of lowest eigenvalue=
2.25733396]
[ 0.46982945
0.74734234
The other eigenvectors could be computed in the same way.
113
0.46982945]
0.8
0.6
relative magnitude
0.4
0.2
0.0
−0.2
−0.4
−0.6
lambda_1
lambda_2
lambda_3
lambda_4
−0.8
1
2
3
4
index
5
6
7
Figure 37: A plot of eigenvectors corresponding to the four lowest eigenvectors in Problem 14.
Problem 14
Note that for this system A is tridiagonal so we can use eigenvals3 with N = 4 to get
the four smallest eigenvalues. We can then use inversePower3 on each to get the eigenvectors corresponding to a given eigenvalue. This problem is worked in the python code
c 9 2 p 14.py. When we run that code we get the following
evalues=
[
4.99017072e-02
7.93335030e+01
1.79505583e+02
5.70789691e+02]
In addition, we generate the plot given in Figure 37. Notice that on that plot each eigenvector
is either “on” (has a nonzero value) for the masses to the left or the right of the mass that
has the small coupling constant k4 . This is because this small constant decouples the masses
at this spring and the system of masses to the left move “independently” from the system
of masses on the right.
Problem 15
Part (a): From the given statement we have an eigenvalue problem of the form
Ay = λBy .
To transform this to a standard form eigenvector problem multiply both sides by the matrix
B −1 to get B −1 Ay = λy. From the given form of B we see that B −1 is a diagonal matrix
with ones on the diagonal and the value of 2 in the (n, n)th location. As B is diagonal so
is B −1 and the product B −1 A is still tridiagonal. Given this tridiagonal matrix we can use
114
eigenvals3 method with N = 1 to find the smallest eigenvalue and then inversePower3
to compute its eigenvector.
We now compute the expression in the book for the value of ω1 . From the equation
ρ
y,
y ′′ = −ω 2
E
we know that there is a solution y(x) of the form
r
ρ
y(x) = A sin ω
x ,
E
with the value of ω determined by the boundary condition y ′(L) = 0. Note that the above
expression for y(x) satisfies the boundary condition of y(0) = 0 already. From the above
form we have that the first derivative of y(x) is given by
r
r
ρ
ρ
′
y (x) = A
ω cos ω
x .
E
E
Thus to have y ′ (L) = 0 we must have
r
ρ
π
L = (2k − 1) ,
ωk
E
2
for k ≥ 1. If we want the smallest eigenvalue then we should take k = 1 and get
s
π 1
E
ω1 =
,
2 L
ρ
(18)
as claimed in the book. Numerically we are actually computing the eigenvalues of the linear
system which we denote as λi . These are related to the eigenvalues of the continuous system
via
2
ωi L ρ .
λi =
n
E
To check our numerical algorithm we will compute the smallest eigenvalue of the linear
system when L = ρ = E = 1. This gives
ω 2
1
λ1 =
,
n
Solving for ω1 we get that
p
ω 1 = n λ1 .
Part (b): This problem is worked in the python code c 9 2 p 15.py. When we run that
code we get the following
n=
10:
n= 100:
n= 500:
n= 750:
n= 1000:
lambda_1=
lambda_1=
lambda_1=
lambda_1=
lambda_1=
0.024623;
0.000247;
0.000010;
0.000004;
0.000002;
omega_1_approx=
omega_1_approx=
omega_1_approx=
omega_1_approx=
omega_1_approx=
115
1.569182;
1.570780;
1.570796;
1.570796;
1.570796;
omega_1_truth=
omega_1_truth=
omega_1_truth=
omega_1_truth=
omega_1_truth=
1.570796
1.570796
1.570796
1.570796
1.570796
Note in the above python code I am using full matrices and numpy routines to compute the
eigenvalues of the given system and then comparing it with the truth. I attempted to do this
with the routines from this book but found stability issues when n got larger. For example
running with compute will full matrices set to False we get
n=
10: lambda_1=
0.024623; omega_1_approx=
1.569182; omega_1_truth=
n= 100: lambda_1=
0.000729; omega_1_approx=
2.700728; omega_1_truth=
n= 500: lambda_1=
0.000066; omega_1_approx=
4.071551; omega_1_truth=
n= 750: lambda_1=
0.000029; omega_1_approx=
4.071618; omega_1_truth=
ridder.py:17: RuntimeWarning: overflow encountered in double_scalars
s = sqrt(fc**2 - fa*fb)
n= 1000: lambda_1=
0.000011; omega_1_approx=
3.324626; omega_1_truth=
1.570796
1.570796
1.570796
1.570796
1.570796
Which shows that we get approximately correct results for small n but the results get worse
as n increase. I didn’t have time to look into what was the cause of this instability (it looks
to be happening in ridder.py).
Problem 16
This problem is worked in the python code c 9 2 p 16.py. When we run that code we get
the following
n=
25; three lowest matrix eigenvalues (lambda)= [ 0.09577011
0.14680479
0.16464815]
Problem 17
This problem is worked in the python code c 9 2 p 17.py. When we run that code we get
the following
evalues= [ 0.00379334 0.01515898 0.0340538
0.06040613 0.094116
eigenvalues computed with np.linalg.eig=
[0.00394654, 0.01458225, 0.0354255, 0.05811637, 0.09788697]
Notice that the results are close to that when I use np.linalg.eig.
116
]
Chapter 10 (Introduction to Optimization)
Problem Set 10.1
Problem 1
This problem is worked in the python code c 10 1 p 1.py. When we run that code we get
the following
x= 0.890898716918
f(x)= -0.25
The analytical solution is when the derivative equals zero. We can show that this equals
1/6
1
≈ 0.8908987181403393 ,
2
in agreement with the number we computed above.
Problem 2
This problem is worked in the python code c 10 1 p 2.py. When we run that code we get
the following
x= 3.53137299707
f(x)= -3.5819473202
Problem 3
This problem is worked in the python code c 10 1 p 3.py. When we run that code we get
the following
x= 1.65215838454
f(x)= -0.844104332014
Problem 4
This problem is worked in the python code c 10 1 p 4.py. When we run that code we get
the following
117
x= 2.28712866065
f(x)= [-672.13578501]
Problem 5
This problem is worked in the python code c 10 1 p 5.py. When we run that code we get
the following
x= 15.9999998973
f(x)= 731.149856985
Problem 6
This problem is worked in the python code c 10 1 p 6.py. When we run that code we get
the following
Constraint Lambda= 1000.0
Minimium point= [ 0.599801
0.4007982]
Constraints satisfied (all negative or zero)?
Function value= 0.519600678364
Constraint Lambda= 10000.0
Minimium point= [ 0.59998001 0.40007998]
Constraints satisfied (all negative or zero)?
Function value= 0.519960006798
Constraint Lambda= 100000.0
Minimium point= [ 0.599998 0.400008]
Constraints satisfied (all negative or zero)?
Function value= 0.519996000068
Constraint Lambda= 1000000.0
Minimium point= [ 0.5999998 0.4000008]
Constraints satisfied (all negative or zero)?
Function value= 0.519999600001
(0.00059920183172623709, 0.00019900308009324075)
(5.9991990270180651e-05, 1.9990085695975601e-05)
(5.9999453041470474e-06, 1.9999091556144322e-06)
(6.0001099733142382e-07, 2.0000280975818185e-07)
Here you see the sequence of larger and larger values of the multiplier λ and the sequence of
solution values that result when we use the output of one solve as the starting location for
the next solve. From this output it looks like the optimal point x is converging to the value
of (0.6, 0.4) which would satisfy both constraints exactly.
Problem 7
This problem is worked in the python code c 10 1 p 7.py. When we run that code we get
the following
118
Constraint Lambda= 1.0
Minimium point= [-0.00231481 0.02777778]
Constraints satisfied (all negative or zero)?
Function value= -1.07167352538e-05
Constraint Lambda= 10.0
Minimium point= [-0.00231481 0.02777778]
Constraints satisfied (all negative or zero)?
Function value= -1.07167352538e-05
(-0.027777778134398575,)
(-0.02777777806009778,)
Notice that the constraints in this problem are not active and thus increasing the value of λ
does not change the solution found.
Problem 8
Notice that the solution to the above problem also satisfies the constraint requirements of
this problem and thus is a solution to this problem also. This problem has a larger search
domain (y ≥ −2 vs. y ≥ 0) and so could have a value of x that results in a smaller function
value. This problem is worked in the python code c 10 1 p 8.py. When we run that code
we get the following
Constraint Lambda= 1.0
Minimium point= [-0.00231481 0.02777778]
Constraints satisfied (all negative or zero)?
Function value= -1.07167352538e-05
Constraint Lambda= 10.0
Minimium point= [-0.00231481 0.02777778]
Constraints satisfied (all negative or zero)?
Function value= -1.07167352538e-05
(-2.0277777781343986,)
(-2.0277777780600976,)
Which is the same solution as found in Problem 8 indicating that increasing the problem
domain did not change the optimal value.
Problem 9
The distance d from any point (x, y) to the point (1, 2) is given by
d2 = (x − 1)2 + (y − 2)2 .
If the point (x, y) must lie on the equation y = x2 then the above can be written as
d2 = (x − 1)2 + (x2 − 2)2 .
119
This is the objective we want to minimize. Note that this is now a one-dimensional problem.
This problem is worked in the python code c 10 1 p 9.py. When we run that code we get
the following
x= 1.36602540313
y=x^2= 1.86602540201
f(x)= 0.151923788647
Indicating that the point on the parabola y = x2 that is closest to (1, 2) is the point
(1.36602540313, 1.86602540201).
Problem 10
Given the figure in the text to determine the y location of the centroid we will break the
total region down into three smaller rectangular regions each of which we can easily compute
the y centroid of. Once we have the centroids of these three rectangles we will determine the
y location of the centroid of the entire object using the formula
P
i Ciy Ai
Cy = P
,
(19)
i Ai
where Ai is the area of the ith region and Ciy is the y centroid of the ith region. The three
regions will be
• the leftmost “tall” rectangle with height 0.4 and width 0.1. This has C1y = 0.2 and
A1 = 0.4(0.1) = 0.04.
• the middle “squashed” rectangle with height 0.4 − x and width 0.2. This has C2y =
1
(0.4 − x) and A2 = 0.2(0.4 − x).
2
• the rightmost “tall” rectangle with height 0.4 and width 0.1. This has C3y = 0.2 and
A3 = 0.4(0.1) = 0.04.
Then with the above decomposition the total area A is
A=
3
X
i=1
and the sum
P3
i=1
3
X
i=1
Ai = 0.08 + 0.08 − 0.2x = 0.16 − 0.2x ,
Ciy Ai is
Ciy Ai = 0.008 + 0.1(0.4 − x)2 + 0.008 = 0.016 + 0.1(0.4 − x)2 .
120
Using these in Equation 19, the y location of the centroid of the entire figure (in terms of x)
is given by
0.016 + 0.1(0.4 − x)2
Cy =
.
0.16 − 0.2x
It is this expression we desire to minimize as a function of x. the python code c 10 1 p 10.py.
When we run that code we get the following
x= 0.234314572512
f(x)= 0.165685424949
Problem 11
If the depth of the water is x then the centroid of the water plus cylindrical vessel is given
by
M(0.43H) + ((πr 2 )xρ)(x/2)
.
M + πr 2 ρx
This is the function we seek to minimize. This problem is worked in the python code
c 10 1 p 11.py. When we run that code we get the following
x= 0.278015691247
f(x)= 0.278015691338
Problem 12
We want to have the volume of the box V fixed at the value 1.0 once folded. In terms of the
variables of the problem this means that
V = ba2 = 1 ,
since b ends up as the height of the box when folded. The area of the cardboard box (once
we remove the four corner squares of size b × b) is given by
(a + 2b)2 − 4b2 .
Thus our problem is find (a, b) that minimize
(a + 2b)2 − 4b2 ,
subject to ba2 = 1. This last constraint means that b =
objective function becomes
2
2
4
a+ 2 − 4 .
a
a
121
1
a2
and using that expression our
Since there is only the variable a in the above we have a one-dimensional problem. In fact,
we can find the minimum of this expression analytically. Taking the a derivative of this
function gives
2
4
16
2 a+ 2
1− 3 + 5 .
a
a
a
We would want to set this expression equal to zero and solve for a. Doing this we get a
polynomial in the variable a given by
a3 − 2 = 0 or a = 21/3 = 1.25992105 .
This problem is worked numerically in the python code c 10 1 p 12.py. When we run that
code we get the following
x= 1.25992104133
f(x)= 4.7622031559
The same as our analytic result above.
Problem 13
For this problem we want to minimize the given expression for V as a function of u and v.
This problem is worked in the python code c 10 1 p 13.py. When we run that code we get
the following
Minimium point (mm)= [ 5.21061757
Function value= -0.106693006162
numIter= 4
28.37508795]
Problem 14
From the diagram it looks like θ = π6 is decent initial guess at the optimal value of θ. To
get a valid initial guess for A we assume θ = π6 and find what value of A would make our
constraints hold with equality. Since there are two constraints we can do this for each of
them. When we do this we find values of Aec (for “equality constraint”) given by
First constraint (without A)= 50000.000000, target= 150000000.000000, A_ec=
0.000333
Second constraint (without A)=
0.000001, target=
0.005000, A_ec=
0.000231
Initial guess for A=
0.000282
Our initial guess for A is taken to be the average of the two “equity constraint” values. Since
our objective function has A in the denominator if A is ever negative the entire expression
122
will becomes negative and our minimization will march off to −∞. To prevent this from
happening we add the constraint A > 0. To make sure that that this constraint is not
violated we need to specify a relatively large value of λ. This problem is worked in the
python code c 10 1 p 14.py. When we run that code we get the following
Minimium point= [ 6.31963854e-01
2.82136567e-04]
Constraints satisfied (all negative or zero)? [-0.0003, -0.632, -0.9388, 0.0, -0.0019]
Function value= 0.00139867355159
Constraint Lambda= 1000000.0
Minimium point= [ 6.31963854e-01
2.82136567e-04]
Constraints satisfied (all negative or zero)? [-0.0003, -0.632, -0.9388, 0.0, -0.0019]
Function value= 0.00139867355159
Thus converting the above angle in radians to degrees we see that the optimal value found
for θ in degrees is 36.20886. Here we see that the constraint on σ is tight.
Problem 15
In this problem we now take δ ≤ 2.5 10−3 meters. This problem is worked in the python
code c 10 1 p 15.py. When we run that code we get the following
Constraint Lambda= 100000.0
Minimium point= [ 7.85398153e-01
2.82482060e-04]
Constraints satisfied (all negative or zero)? [-0.0003, -0.7854, -0.7854, -24840407.7816, 0.0]
Function value= 0.00159897861142
Constraint Lambda= 1000000.0
Minimium point= [ 7.85398164e-01
2.82806522e-04]
Constraints satisfied (all negative or zero)? [-0.0003, -0.7854, -0.7854, -24984004.0909, 0.0]
Function value= 0.0015998976262
Here we see that the constraint on σ is no longer tight but the constraint on δ is now tight.
Problem 16
For this problem we added the constraints that r1 > 0 and r2 > 0 and need a large value of
λ to make sure that these constraints are satisfied. This problem is worked in the python
code c 10 1 p 16.py. When we run that code we get the following
Constraint Lambda= 1000.0
Minimium point (mm)= [ 53.69286082 41.35669939]
Constraints satisfied (all negative or zero)? [-0.0537, -0.0414, -15490689.7861, -0.0, 0.0001]
Function value= 0.0144463267945
Constraint Lambda= 10000.0
Minimium point (mm)= [ 53.77874435 41.35669939]
Constraints satisfied (all negative or zero)? [-0.0538, -0.0414, -16277585.6, -0.0, 0.0]
Function value= 0.0144609102567
123
From the above output we see that the first stress constraint (the one on σ1 ) is not tight but
the one on σ2 and δ is.
Problem 17
To find the analytic solution we set the derivatives of F (x, y, z) equal to zero as
∂F
= 4x + y + z = 0
∂x
∂F
= 6y + x − 2 = 0
∂y
∂F
= 2z + x = 0 .
∂z
As a matrix system this is

   
4 1 1
x
0
 1 6 0  y  =  2  .
1 0 2
z
0
We can solve this linear system for x, y, z to get the analytic solution for the minimum. We
can also use a routine like powell to compute a numerical solution. This problem is worked
in the python code c 10 1 p 17.py. When we run that code we get the following
analytic solution x= [[-0.1
Minimium point= [-0.1
0.35
Function value= -0.35
0.35 0.05]]
0.05]
Problem 18
For this problem we again have to add constraints that the dimensions are positive and
specify a large value of λ so that these constraints are enforced. This problem is worked in
the python code c 10 1 p 18.py. When we run that code we get the following
Constraint Lambda= 1000.0
Minimium point= [ 0.752725
0.33662887 0.67325772]
Constraints satisfied (all negative or zero)? [-0.7527, -0.3366, -0.6733, 0.0013]
Function value= 3.981989384
Constraint Lambda= 10000.0
Minimium point= [ 0.75303353 0.33673959 0.67356016]
Constraints satisfied (all negative or zero)? [-0.753, -0.3367, -0.6736, 0.0001]
Function value= 3.98357718769
124
Problem 19
This problem is worked in the python code c 10 1 p 19.py. When we run that code we get
the following
Constraint Lambda= 1.0
Minimium point= [ 0.59649121 1.
0.75521228]
Constraints satisfied (all negative or zero)? [-8546.4659, -5469.7381, 0.0001]
Function value= [ 2.35170351]
Constraint Lambda= 10.0
Minimium point= [ 0.53333333 1.
0.73333333]
Constraints satisfied (all negative or zero)? [0.0, 0.0, 0.0]
Function value= [ 2.26666667]
Problem 20
This problem is worked in the python code c 10 1 p 20.py. When we run that code we get
the following
Constraint Lambda= 1000000.0
Minimium point (degrees)= [ 22.91215684
1.53852011 -29.54687013]
Constraints satisfied (all negative or zero)? [0.0253, 0.0143]
Function value= -23723.9572909
Constraint Lambda= 10000000.0
Minimium point (degrees)= [ 21.52804395
1.27824854 -28.18821937]
Constraints satisfied (all negative or zero)? [0.0027, 0.0014]
Function value= -22928.3090388
Constraint Lambda= 100000000.0
Minimium point (degrees)= [ 21.37733298
1.251816
-28.03662506]
Constraints satisfied (all negative or zero)? [0.0003, 0.0001]
Function value= -22844.203737
125
References
[1] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning. Springer,
New York, 2001.
126