Math 466 Numerical Analysis

Math 466 Numerical Analysis
Time and place:

T, 8:30-9:20, C-401 Padelford
Instructor:
Jim Morrow
o ce hours:
MW 8:30-9:20 C439 Padelford
Prerequisites:
Math 465 and Math 307
Textbooks:
Numerical Analysis by Johnson and Riess
This course will introduce some numerical techniques for solving di erential equations.
It will be expected that you spend some time reading the textbook. You can tell from the
homework assignments which sections you should be reading.
The homework assignments are as follows:
due date
assignment
April 6
x7.1: 8,11; x7.2.2: 4,10,13,14
April 13
x7.2.3: 5,6,8,9;
April 20
x7.2.5: 2,4,6,7
April 27
x7.3: 1,3,7
May 4
x7.4: 2,7,8,9; 7.5.1: 5,13
May 11
x7.5.3: 1c,2c,4,10; 7.5.4: 3,4,5
May 18
x7.6.3: 6,7,8
May 25
x8.1: 5,8
June 1
x8.2.1: 3,4
Project # 1 (due April 20)
Write a program RK4 which uses the fourth order Runge-Kutta method (Johnson &
Riess, p369-371) to calculate a numerical solution to an initial value problem:
y0 = f (x; y); y(a) = ya

RK4 should take as inputs the following arguments:
initial x value a
nal x value b
initial y value ya
positive integer m, the number of subdivisions
The step-size is then h = (b a)=m. There is also a parameter VALMAX, which may
be taken to be 1.e4. RK4 should call on a subroutine for the function evaluations. This
subroutine func(x; y; z) takes the values x; y and places the computed value in z.
RK4 should return an error message if some computed y or y0 has an absolute value
that exceeds VALMAX. The program should return the computed values (x ; y ) at m + 1
equally spaced points in a; b], or as many points as are computed before VALMAX is
exceeded. Have the results sent to an external le.
i
1. Run the program with f (x; y) = x 100y, and a = 0; b = 1; ya = 1. Let m = 2 , where

k = 2; 3; : : : ; 8. Compare the computed solution values with the values obtained by solving
k
the di erential equation exactly. Plot the computed solutions and the true solution on a
single graph.
2. Let f (x; y) = y(2 y). Run the program with a = 0; b = 4; k = 8, for the three initial
y values ya = :1; 1; 3, and let m = 28 . Plot the results on a single graph.
3. Let f (x; y) = x2 + y2 . Let a = 0; b = 1; ya = 1, and let m = 2 , where k = 2; : : : ; 8.
Plot the solutions and explain the results.
k
Project # 2 (due May 11)

Write a program that uses the Heun method to solve a system of ODE's:
u0 = F(x; u)
where u and F are n-vectors. Take the suggestion of Johnson & Riess, p396, and write a
subroutine func(x; u; y), which puts the value of F(x; u) in y. The program should take
as inputs the following arguments:
integer n 10, where n is the dimension of the system
initial x value
nal x value
initial values u1; u2; : : : ; u
positive integer m, where the number of subdivisions is m
The step-size is then h = (b a)=m. There is also a parameter VALMAX, which may
be taken to be 1.e6.
The program should place the values u(x ) = (u1(x ); u2 (x ); : : : ; u (x )), calculated
at the intermediate points x , in an external le or print an error message if some calculated
value of ku(x )k or ky(x )k exceeds VALMAX.
n
1. Run your program on the system:
u01 = 30u1 + 28u2

2u2
u02 =
with 0 x 4, u1(0) = 2; u2(0) = 1 and k = 3; : : : ; 8. Plot the result for m = 28 in the
(u1; u2)-plane. Do the same thing for:
u01 = 30u1 + 28u2

u02 =
2u2
which values of m give reasonable answers? 2. Convert the second order ODE
y00 + xy0 + x2 y = 0
to a rst order system, and then run your program on the system to nd the approximate
value of y(4), given the initial conditions y(0) = 1; y0 (0) = 0. Let m = 2 where k =
3; : : : ; 8. Plot the result for k = 8 in the (x; y)-plane.
k
Project # 3 (due May 18)

Use three-step Adams-Bashforth as a predictor and two-step Adams-Moulton (formulas 7.47a and 7.47b of Johnson & Riess) to write a program to solve y0 = f (x; y). Use your
RK4 routine to compute the starting values. The program should take as inputs:
initial x value xa
nal x value xb
initial y value ya
number of subdivisions m
There should be a parameter VALMAX set equal to 1.e5.
The program should return the computed values (x ; y ) or an error message if jy j or
jf (x ; y )j exceeds VALMAX. Use this program as follows:
i
1. Let f (x; y) = y=2; xa = 0; xb = 2; ya = 1; m = 2 ; k = 3; : : : ; 8. Plot the results.

2. Let f (x; y) = 80(sin(x) y); xa = 0; xb = 3; ya = 0; m = 2 ; k = 3; 5; 7; 9. Plot
the results. Explain the results. Find experimentally the smallest m that can be used to
give a reasonable approximation to the true solution. Solve this equation using the RK4
method of project # 1.
3. The global error is usually not the sum of the local errors. Take the initial value
problem of part 1 and k = 8, and compute the di erence between the sum of the local
errors and the global error.
4. The relative global error should be closer to the sum of the relative local errors.
Use part 1 again and compute the di erence between the relative global error and the sum
of the relative local errors.
Project # 4 (Extra Credit, due June 1)
k
1. A disposal plant has an ejector which satis es the ODE:

80y00 + 20(y0 )2 + y = 60; y(0) = 0; y0 (0) = 0
where y(t) is the height of the uid in the ejector chamber at time t. Convert this equation to a rst order system of ODE's. Using stepsize h = 2 6, compute a numerical
approximation, for 0 t 4. Plot the results.
2. (Lorenz's example) Consider the system
X 0 = 10X + 10Y
Y 0 = XZ + rX Y
Z 0 = XY 8Z=3
with initial conditions X (0) = 0; Y (0) = 1; Z (0) = 0. Solve this system on 0 t 30,
using the method of your choice with a stepsize of :01, with r assuming the values 5; 20; 28.
In each case make the following plots: (t; Y ); (X; Y ); (Y; Z ).
3. Solve the system
s0 = cd
c0 = ds
d0 = sc=3
with initial conditions s(0) = 0; c(0) = d(0) = 1. This is the system satis ed by the
angular momentum of a rigid body with moments of inertia (1; 2; 3) (an asymmetric top).
The solution is periodic in t. Find the period P by nding the rst P > 0 for which
s(P ) = 0; c(P ) = d(P ) = 1. This period is also given by an integral
4
Z1
dx
(1 x2)(1 x2 =3)
Use a Gaussian quadrature routine to compute this integral. There is yet another way to
compute P (also due to Gauss) using the arithmetic- p
geometric mean. Let a0 = 1; b0 =
p
2=3. For n 0, let a +1 = (a + b )=2; b +1 = a b . The sequences fa g; fb g
converge very rapidly to a common limit L. Gauss showed that P = 2 =L.
Project # 5 (Extra Credit, due June 1)
n
n n
1. Use nite di erences to approximate the solution of Laplace's equation u + u = 0

in the region R = f(x; y) : 0 x 1; 0 y 1g, with the boundary conditions u(x; 0) =
0; u(0; y) = 0; u(x; 1) = 10 sin x; u(1; y) = 0.
(a) Use the Gauss-Seidel method to solve the matrix equation that arises. Continue the
iteration until two successive approximations di er by less than TOL or until MAX iterations have occurred. Take TOL = :01; MAX = 200. Let the step size h = =N , where
N takes on the values 10; 20; 50. For this problem nd the number of iterations required
for each step size. Do not print out the solution.
(b) Do the this problem (with N = 50) using SOR with the acceleration parameter !
assuming the values 1:0; 1:1; : : : ; 1:9 as well as the optimal value ! = 2=(1 + sin( =50)).
In each case nd the number of iterations until TOL is obtained.
xx
yy
opt
2. Consider the heat equation u = u ; 0 x 1; t 0, with boundary conditions

u(x; 0) = f (x); u(0; t) = 0 = u(1; t).
(a) Solve the equation by the direct method. Let h = 1=M be the step size in the x
direction. Let the t step size be k = rh2 . Let the nal t value be c. The program
should return the values u(x ; c); j = 0; : : : ; M . Run the program with the data: f (x) =
sin x + 4 sin 4 x; M = 20, and c = :0625; :125; :25, and let k = 1=400 (r = 1), and
k = 1=1200 (r = 1=3). Plot the results.
(b) Solve the equation with the implicit method where k = 1=40; 1=100; 1=400 and c =
:0625; :125; :25. Plot the results for k = 1=40.
(c) Solve the equation with the Crank-Nicholson method, taking k = 1=40; 1=100 and
c = :0625; :125; :25. Plot the results for k = 1=40.
t
xx

Math 466 Numerical Analysis

Uploaded by

Copyright:

Available Formats

Math 466 Numerical Analysis

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Math 466 Numerical Analysis

Uploaded by

Copyright:

Available Formats

Math 466 Numerical Analysis

Time and place:

y0 = f (x; y); y(a) = ya

1. Run the program with f (x; y) = x 100y, and a = 0; b = 1; ya = 1. Let m = 2 , where

Project # 2 (due May 11)

1. Run your program on the system:

u01 = 30u1 + 28u2

u01 = 30u1 + 28u2

Project # 3 (due May 18)

1. Let f (x; y) = y=2; xa = 0; xb = 2; ya = 1; m = 2 ; k = 3; : : : ; 8. Plot the results.

1. A disposal plant has an ejector which satis es the ODE:

3. Solve the system

1. Use nite di erences to approximate the solution of Laplace's equation u + u = 0

2. Consider the heat equation u = u ; 0 x 1; t 0, with boundary conditions

You might also like