There are 18 questions in full, however, only FOUR questions are for your coursework, which will be chosen (determined) by the computer (by the method of a random choice). There is no mark if you do the questions selected by yourself.
1. (a) Avoiding large errors, overflow and underflow. For 100 values of x over the interval [1014
, 1016], evaluate the following expression
y =
p
2x
2 + 1 − 1
and then develop a new expression which is mathematically equivalent, plot them, and based on the graphs, tell which is better in terms of resisting the loss of significance.
(b) Implement Newton’s method. Write m-file with a function like function x=run_Newton(f, df, x0, N, error) where f is the name of the function, and df is its derivative. Run your code to find zeroes fo the following functions
(a) f(x) = tan x − x
(b) f(x) = (x − 1)7
(c) f(x) = x
7 − 7x
6 + 21x
5 − 35x
4 + 35x
3 − 21x
2 + 7x − 1
(d) f(x) = x
2 + 1
(Note: the last function has no zeroes, please give your discussions)
2. Write Matlab codes to find the roots of the functions.
(a) The volume V of liquid in a spherical tank of radius r is related to the depth h of the liquid by
V =
πh2
(3r − h)
3
i. Determine h by using Newton’s method given r = 1m and V = 0.75m3.
ii. For such a spherical tank, it is possible to develop the following two fixedpoint formulas h =
r
h3 + (3V /π)
3r
and
h =
3
s
3
rh2 −
V
π
If r = 1m and V = 0.5m3
, determine whether either of these is stable,
and the range of initial guesses for which they are stable. Please give your discussions.
(b) Consider the function
f(x) = 14xex−2 − 12e
x−2 − 7x
3 + 20x
2 − 26x + 12
on the interval [0, 3]
2
i. Plot the function on the interval
ii. Apply Newton’s method to find both roots
iii. Determine for which roots Newton’s method convergences linearly and for which the convergence is quadratic
3. Find the roots for the following problems.
(a) Here is a cubic polynomial with three closely spaced real roots:
p(x) = 816x
3 − 3835x
2 + 6000x − 3125 = 0
i. What are the exact roots of p(x) = 0?
ii. Plot p(x) for 1.43 ≤ x ≤ 1.71. Show the location of the three roots.
iii. Starting with x0 = 1.5, what does Newton’s method do?
iv. Starting with x0 = 1 and x1 = 2, what does bisection do?
Please give your discussion for these questions.
(b) The idea gas law for a gas at low temperature and pressure is P V = nRT where P is pressure (in atm), V is volume (in L), T is temperature (in K), n is the number of moles of the gas, and R = 0.0820578 is the molar gas constant.
Van der Waal’s equation
P +
n
2a
V 2
(V − nb) = nRT covers the non-ideal case where these assumptions do not hold. Use the ideal gas law to compute an initial guess, followed by Newton’s method applied to van der Waals’ equation to find the volume of one mole of oxygen at 320 deg.K. and a pressure of 15 atm. For oxygen, a = 1.36 L 2 atm/mole2 and b = 0.003183 L/mole. State your initial guess and solution with three significant digits.
(c) Use the same code, to find the volume of 1 mole of benzene vapor at 700 deg.K. under a pressure of 20 atm. For benzene, a = 18.0L
2 .atm/mole2 and b = 0.1154
L/mole
4. Numerical integrations
(a) Apply the Simpson’s method to the integrals with 20 < n < 200
(a)
Z 1
0
e
x
2
dx, (b)
Z π
2
0
e
x − 1
sin x
dx (c)
Z π
0
e
cos x
dx
You must choose at least three different ns to calculate the solutions and make your discussions for your numerical results.
(b) The composite trapezoid rule with n equally spaced points is Tn(f) = h
2
f(a) + h
nX−2
k=1
f(a + kh) + h
2
f(b)
where
h =
b − a
n − 1
3
Use Tn(f) with various values of n to compute π by approximating
π =
Z 1
−1
2
1 + x
2
dx
How does the accuracy vary with n. Please give your discussion.
5. Solve linear systems
Write Matlab codes for both Jacobi and Gauss-Seidel methods, respectively, to solve the following equations. Make sure at the beginning you solve each equation for the variable that has the largest coefficient, why?
(a) Apply your codes to the system of equations
x1 + 9x2 − 2x3 = 36
2x1 − x2 + 8x3 = 121
6x1 + x2 + x3 = 107 with an initial
x1 = 1, x2 = 1, x3 = 1
(b) Apply your codes to solve the system
A(t)x = b
starting from x = [0, 0, 0]
T, where
A(t) =
?
?
1 t t
t 1 t
t t 1
?
? , b =
?
?
2
2
2
?
?
For t = 0.2, 0.5, 0.8, 0.9, 0.99, determine the number of steps (or iterations) to obtain the solution to 6 significant digits in the computations. Plot the number of steps as functions of t and make your comments.
(c) Solve the sparse system within six correct decimal places for n = 100, n = 1000 and n = 10000 by using your codes which you made for solving (a) and (b). The correct solution is xe = [1, 1, · · · , 1]T . Report the number of step iterations and the error max |xn − xe| for all n, here xn is the numerical solutions and xe is the exact solution. The system is
?
??????
3 −1
−1 3 −1
.
.
.
.
.
.
.
.
.
−1 3 −1
−1 3
?
??????
?
???????
x1
x2
.
.
.
.
.
.
xn
?
???????
=
?
??????
2
1
.
.
.
1
2
?
??????
6. Linear algebraic equations can arise in the solution of difference equations. For example, the following differential equation derives from a heat balance for a long, thin rod
d
2T
dx2
+ κ(Ta − T) = 0
4
where T is temperature (0C), x is the distance along the rod (m), and κ is a heat transfer coefficient between the rod and the ambient air (m−2 ), and Ta is the temperature of the surrounding air. This equation can be transformed into a set of linear algebraic equations by using a finite difference approximation for the second derivative d
2T
dx2
=
Ti+1 − 2Ti + Ti−1
h2
where Ti designates the temperature at node xi
. This approximation can be written
as
−Ti−1 + (2 + κh2
)Ti − Ti+1 = κh2Ta
(a) Write this equation for each of the interior nodes of the rod resulting in a tridiagonal system of equations. The first and the last nodes at the rod ends are fixed by the boundary conditions.
(b) Develop an analytical solution for the equation with 10m rod with Ta = 20,,
T(x = 0) = 40, T(x = 10) = 200 and κ = 0.02.
(c) Develop numerical solutions for the same parameters values for different h = 2, h = 1, h = 0.1, h = 0.01 and h = 0.0001. Plot these results and compare with the analytical solution and give your discussions.
7. A cylindrical pipe has a hot fluid flowing through it. Because the pressure is very high, the walls of the pipe are thick. For such a situation, the temperature distribution in the metal wall related to radial distance can be modelled by a second order differential equation
r
d
2
θ
dr2
+
dθ
dr = 0
where θ is temperature, and r is the radial distance from the centreline. In the following questions, you must use both the shooting method and finite difference method to solve the above equation, and compare the solutions with different mesh steps, h, by each method. Plot solutions and make your judgment on solutions for each method.
(a) Solve for the temperature θ(r) within a pipe whose radius is 1 cm and whose outer radius is 2 cm if the fluid is at 500oC and the temperature of the outer circumference is 25oC.
(b) The pipe is insulated to reduce the heat loss. The insulation used has properties such that the gradient dθ dr at the outer circumference is proportional to the difference in temperature from the outer wall to the surroundings, which can be modelled by
dθ
dr
r=2
= 0.083(θ(2) − 20)
Solve the equation with this boundary condition.
8. The governing equation for a projectile shot vertically upward is
m
d
2y
dt2
= −mg − C|V |V, y(0) = 0, y′
(0) = V0
where m is the mass of the projectile (kg), y(t) is the height (m), g is the acceleration
of gravity (9.80665m/s2
), C is an aerodynamic drag parameter, and V = dy/dt is
5
the velocity. For m = 10kg, C = 0.1Ns2/m2
, and V0 = 500.0m/s, calculate (a)
the maximum height attained by the projectile, (b) the time required to reach the maximum height, and (c) the time required to return to the original elevation.
9. A cable hung at its two ends as shown in Fig. 1 by its own weight will have a catenary shape described by the equation
y =
Tx
w
cosh wx
Tx
− 1
(1)
where w is the weight per unit length and Tx is the horizontal, x-component of the tension of the cable. Equation (1) is for the case when both w and Tx are constant throughout the cable. In fact, Equation (1) is the solution of the differential equation
d
2y
dx2
=
w
Tx
ds
dx =
w
Tx
"
1 +
dy
dx2
#1/2
(2)
where s is a variable along the length of the cable.
(a) Applying the Runge-Kutta method to solve the equation (2), where w = 0.12 kN/m, xA = yA = 0, xB = 200 m, and yB = 50 m. Let the initial conditions be y = 0 at x = xA = 0, iterate Tx value until yB is within 99.9% of 50 m.
(b) How could the problem be solved if xA = 100 m and yA = 25 m by application of the Runge-Kutta method (w remains equal to 0.12 kN/m)?
(c) Making a discussion about numerical results with different mesh steps and with other methods if possible.
10. The equation
y
′ = 1 + y
2
, y(0) = 0
has the analytical solution y(t) = tan t. The tangent function is infinite at t = π/2.
(a) Write your program by using both Euler and Runge-Kutta method to solve this initial value problem between t = 0 and t = 1.6, and compare the results of the
program with the analytical solution
(b) What is the behaviour of the Runge-Kutta fourth order method when used between t = 0 and t = 1.6
6 (c) Compare the results from your codes and make comments on the behaviour change with step sizes.
11. Given the following non-linear boundary value problem
y
′′ − 18y
2 = 0, y(1) = 1
3
, y(2) = 1
12
(a) Use the shooting method to approximate solution
(b) Use finite difference to approximate solution
(c) Plot the approximate solutions together with the exact solution y(t) = 1
3t
2 and
discuss your results with both methods
12. Here is a typical steady-state heat flow problem. Consider a thin steel plate to be a 10 × 20 (cm)2
rectangle. If one side of the 10 cm edge is held at 1000C and the other three edges are held at 00C, what are the steady-state temperature at interior points? We can state the problem mathematically in this way if we assume that heat flows only in the x and y directions:
Find u(x, y) (temperature) such that
∂
2u
∂x2
+
∂
2u
∂y2
= 0 (3)
with boundary conditions
u(x, 0) = 0
u(x, 10) = 0
u(0, y) = 0
u(20, y) = 100
We replace the differential equation by a difference equation
1
h2
[ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j ] = 0 (4)
which relates the temperature at the point (xi
, yj ) to the temperature at four neighbouring points, each the distance h away from (xi
, yj ). An approximation of Equation
(3) results when we select a set of such points (these are often called as nodes) and
find the solution to the set of difference equations that result.
(a) If we choose h = 5 cm , find the temperature at interior points.
(b) Write a program to calculate the temperature distribution on interior points with h = 2.5, h = 0.25, h = 0.025 and h = 0.0025 cm. Discuss your solutions and examine the effect of grid size h.
(c) Modified the difference equation (4) so that it permits to solve the equation
∂
2u
∂x2
+
∂
2u
∂y2
= xy(x − 2)(y − 2)
on the region
0 ≤ x ≤ 2, 0 ≤ y ≤ 2
with boundary condition u = 0 on all boundaries except for y = 0, where u = 1.0. Write and run the program with different grid sizes h and discuss your numerical results.
13. If a anti level beam of length L, which bends due to a uniform load of w lb/ft, is also subject to an axial force P at its free end (see Figure 2). The equation of its elastic curve is
EI d
2y
dx2
= P y −
wx2
2
(5)
For this equation, the origin O has been taken at the free end. I is the moment of inertia; here I =
bh3
12
. At the point x = L,
dy
dx = 0 and at x = 0, y = 0.
P
O
X
Y
L
Figure 2: Diagram for a cantilever beam.
(a) Solve the boundary value problem by any numerical method for a wooden beam, 2 in × 4 in × 10 ft, where E = 12 × 105 lb/in2
. Find y versus x when the beam has the 4 in dimension vertical with w = 25 lb/ft and a tension force of P = 500 lb.
(b) And also solve for the deflections if the beam is turned so that the 4 in dimension is horizontal.
(c) If the factor of the radius of curvature has been taken into account, Equation (5) becomes
EI
[1 + (y
′)
2]
3/2
d
2y
dx2
= P y −
wx2
2
(6)
If the defection of the beam is small, the difference is negligible, but in some cases this is not true. Furthermore, if there is considerable bending of the beam, the horizontal distance from the origin to the wall is less than L, the original length of the beam. Solve Equation (6) with the same conditions, and determine by how much the deflections differ from those previous calculations.
14. Consider a laterally insulated metal bar of length L = 1, which satisfies the heat equation
∂u
∂t =
∂
2u
∂x2
Suppose that the ends of the bar are kept at temperature u = 0oC and the temperature in the bar at t = 0, u(0, x) = sin πx.
(a) Using an explicit method with different meshes h to solve the equation and find the numerical solutions, then plot the solutions for t = 0, t = 1, t = 2, t = 3, t = 4 and t = 5 in the same figure.
(b) Using an implicit method with different meshes h to solve the equation and find the numerical solutions, then plot the solutions for t = 0, t = 1, t = 2, t = 3, t = 4 and t = 5 in the same figure.
(c) Using Crank-Nicolson method with different meshes h to solve the equation and find the numerical solutions, then plot the solutions for t = 0, t = 1, t = 2, t = 3, t = 4 and t = 5 in the same figure.
(d) Making your discussion for these three methods with different size of h The exact solution for this problem is
u(t, x) = sin(πx)e
−π
2
t
which can be used to compare with your numerical solutions.
15. Condon and Odishaw (1967) discuss Duffing equation for the flux φ in a transformer. This nonlinear differential equation is
d
2φ
dt2
+ ω
2
0φ + bφ3 =
ω
N
E cos(ωt), φ(0) = φ
(0) = 0
In this equation , E cos ωt is the sinusoidal source voltage and N is the number of turns in the primary winding, where ω0 and b are parameters of the transformer design. Make a plot of φ versus t (and compare to the source voltage) if E = 165, ω = 120π,
N = 600, ω
2
0 = 83 and b = 0.14. For approximate calculations, the nonlinear term bφ3
is sometimes neglected. Evaluate your result to determine whether this makes a significant error in the results.
(Condon E. and Odishaw H., (1967) Handbook of Physics, New York, McGraw-Hill)
16. An 80kg paratrooper is dropped from an airplane at a height of 600m. After 5 second the chute opens. The paratrooper’s height as a function of time, y(t), is given by
d
2y
dt2
= −g +
α(t)
m
y(0) = 600 (m)
y
′
(0) = 0 (m/s)
where g = 9.81m/s2 is the acceleration due to gravity and m = 80kg is the paratrooper’s mass. The air resistance α(t) is proportional to the square of the velocity, with different proportionality constants before and after the chute opens,
α(t) =
K1y
′
(t)
2
, t < 5s
K2y
′
(t)
2
, t ≥ 5s
(a) Find the analytical solution for free fall, K1 = 0 and K2 = 0. At what height does the chute open? How long does it take to reach the ground? What is the impact velocity? Plot the height versus time and label the plot appropriately.
(b) Consider the case K1 =
1
15 . and K2 =
4
15 . At what height does the chute open? How long does it take to reach the ground? What is the impact velocity? Make a plot of the height versus time and label the plot appropriately.
9
17. Apply the shooting method to the boundary value problem
y
′
1 =
4 − y2
t
3
y
′
2 = −e
y1
y1(1) = 0, y2(2) = 0, 1 ≤ t ≤ 2
Note that if the initial condition y2(1) were given, this would be an initial value problem. You may apply the shooting method to determine the unknown y2(1), using matlab routine ODE45 to solve the initial problem. Please give the details about using the shooting method to solve this problem.
The exact solutions are y1(t) = ln t and y2(t) = 2 − t 2/2. You can use it to compare with your numerical solutions.
18. The shooting method for ODE:
(a) Consider the second order differential equation
y
′′ + 2y
′ − 3y = 0, y(0) = e
2
, y(2) = 1
First, find the exact solution for this equation, and then write a matlab code to use the shooting method to find the approximate solutions. In this question, you should give the detail about choosing the initial guesses. Please read lecture note p68-69. You should compare your numerical results with the exact solutions.
(b) Then use your code to solve the following nonlinear second order differential equation
y
′′ + 0.1(y
′
)
2 + 0.6y = sin t, y(0) = 1, y(2) = 2
Clearly, you could not find the exact solution, however, you should plot your
numerical solutions and make a good discussion with different guesses.There are 18 questions in full, however, only FOUR questions are for your coursework,
which will be chosen (determined) by the computer (by the method of a random choice).
There is no mark if you do the questions selected by yourself.
1. (a) Avoiding large errors, overflow and underflow. For 100 values of x over the
interval [1014
, 1016], evaluate the following expression
y =
p
2x
2 + 1 − 1
and then develop a new expression which is mathematically equivalent, plot them, and based on the graphs, tell which is better in terms of resisting the loss of significance.
(b) Implement Newton’s method. Write m-file with a function like function x=run_Newton(f, df, x0, N, error)
where f is the name of the function, and df is its derivative. Run your code to find zeroes fo the following functions
(a) f(x) = tan x − x
(b) f(x) = (x − 1)7
(c) f(x) = x
7 − 7x
6 + 21x
5 − 35x
4 + 35x
3 − 21x
2 + 7x − 1
(d) f(x) = x
2 + 1
(Note: the last function has no zeroes, please give your discussions)
2. Write Matlab codes to find the roots of the functions.
(a) The volume V of liquid in a spherical tank of radius r is related to the depth h
of the liquid by
V =
πh2
(3r − h)
3
i. Determine h by using Newton’s method given r = 1m and V = 0.75m3
ii. For such a spherical tank, it is possible to develop the following two fixedpoint formulas
h =
r
h3 + (3V /π)
3r
and
h =
3
s
3
rh2 −
V
π
If r = 1m and V = 0.5m3
, determine whether either of these is stable, and the range of initial guesses for which they are stable. Please give your discussions.
(b) Consider the function
f(x) = 14xex−2 − 12e
x−2 − 7x
3 + 20x
2 − 26x + 12
on the interval [0, 3]
2
i. Plot the function on the interval
ii. Apply Newton’s method to find both roots
iii. Determine for which roots Newton’s method convergences linearly and for which the convergence is quadratic
3. Find the roots for the following problems.
(a) Here is a cubic polynomial with three closely spaced real roots:
p(x) = 816x
3 − 3835x
2 + 6000x − 3125 = 0
i. What are the exact roots of p(x) = 0?
ii. Plot p(x) for 1.43 ≤ x ≤ 1.71. Show the location of the three roots.
iii. Starting with x0 = 1.5, what does Newton’s method do?
iv. Starting with x0 = 1 and x1 = 2, what does bisection do?
Please give your discussion for these questions.
(b) The idea gas law for a gas at low temperature and pressure is P V = nRT where P is pressure (in atm), V is volume (in L), T is temperature (in K), n is the number of moles of the gas, and R = 0.0820578 is the molar gas constant. Van der Waal’s equation
P +
n
2a
V 2
(V − nb) = nRT covers the non-ideal case where these assumptions do not hold. Use the ideal gas law to compute an initial guess, followed by Newton’s method applied to van der Waals’ equation to find the volume of one mole of oxygen at 320 deg.K. and a pressure of 15 atm. For oxygen, a = 1.36 L 2 atm/mole2 and b = 0.003183 L/mole. State your initial guess and solution with three significant digits.
c) Use the same code, to find the volume of 1 mole of benzene vapor at 700 deg.K. under a pressure of 20 atm. For benzene, a = 18.0L
2
.atm/mole2 and b = 0.1154
L/mole
4. Numerical integrations
(a) Apply the Simpson’s method to the integrals with 20 < n < 200
(a)
Z 1
0
e
x
2
dx, (b)
Z π
2
0
e
x − 1
sin x
dx (c)
Z π
0
e
cos x
dx
You must choose at least three different ns to calculate the solutions and make your discussions for your numerical results.
(b) The composite trapezoid rule with n equally spaced points is
Tn(f) = h
2
f(a) + h
nX−2
k=1
f(a + kh) + h
2
f(b)
where
h =
b − a
n − 1
3
Use Tn(f) with various values of n to compute π by approximating
π =
Z 1
−1
2
1 + x
2
dx
How does the accuracy vary with n. Please give your discussion.
5. Solve linear systems
Write Matlab codes for both Jacobi and Gauss-Seidel methods, respectively, to solve the following equations. Make sure at the beginning you solve each equation for the variable that has the largest coefficient, why?
(a) Apply your codes to the system of equations
x1 + 9x2 − 2x3 = 36
2x1 − x2 + 8x3 = 121
6x1 + x2 + x3 = 107
with an initial
x1 = 1, x2 = 1, x3 = 1
(b) Apply your codes to solve the system
A(t)x = b
starting from x = [0, 0, 0]
T, where
A(t) =
?
?
1 t t
t 1 t
t t 1
?
? , b =
?
?
2
2
2
?
?
For t = 0.2, 0.5, 0.8, 0.9, 0.99, determine the number of steps (or iterations) to obtain the solution to 6 significant digits in the computations. Plot the number of steps as functions of t and make your comments.
(c) Solve the sparse system within six correct decimal places for n = 100, n = 1000 and n = 10000 by using your codes which you made for solving (a) and (b). The correct solution is xe = [1, 1, · · · , 1]T
. Report the number of step iterations and the error max |xn − xe| for all n, here xn is the numerical solutions and xe is the exact solution. The system is
?
??????
3 −1
−1 3 −1
.
.
.
.
.
.
.
.
.
−1 3 −1
−1 3
?
??????
?
???????
x1
x2
.
.
.
.
.
.
xn
?
???????
=
?
??????
2
1
.
.
.
1
2
?
??????
6. Linear algebraic equations can arise in the solution of difference equations. For example, the following differential equation derives from a heat balance for a long, thin rod
d
2T
dx2
+ κ(Ta − T) = 0
4
where T is temperature (0C), x is the distance along the rod (m), and κ is a heat transfer coefficient between the rod and the ambient air (m−2 ), and Ta is the temperature of the surrounding air. This equation can be transformed into a set of linear algebraic equations by using a finite difference approximation for the second derivative
d
2T
dx2
=
Ti+1 − 2Ti + Ti−1
h2
where Ti designates the temperature at node xi . This approximation can be written
as
−Ti−1 + (2 + κh2
)Ti − Ti+1 = κh2Ta
(a) Write this equation for each of the interior nodes of the rod resulting in a tridiagonal system of equations. The first and the last nodes at the rod ends are fixed by the boundary conditions.
(b) Develop an analytical solution for the equation with 10m rod with Ta = 20,, T(x = 0) = 40, T(x = 10) = 200 and κ = 0.02.
(c) Develop numerical solutions for the same parameters values for different h = 2, h = 1, h = 0.1, h = 0.01 and h = 0.0001. Plot these results and compare with the analytical solution and give your discussions.
7. A cylindrical pipe has a hot fluid flowing through it. Because the pressure is very high, the walls of the pipe are thick. For such a situation, the temperature distribution in the metal wall related to radial distance can be modelled by a second order differential equation
r
d
2
θ
dr2
+
dθ
dr = 0
where θ is temperature, and r is the radial distance from the centreline.
In the following questions, you must use both the shooting method and finite difference method to solve the above equation, and compare the solutions with different mesh steps, h, by each method. Plot solutions and make your judgment on solutions for each method.
(a) Solve for the temperature θ(r) within a pipe whose radius is 1 cm and whose outer radius is 2 cm if the fluid is at 500oC and the temperature of the outer circumference is 25oC.
(b) The pipe is insulated to reduce the heat loss. The insulation used has properties such that the gradient dθ dr at the outer circumference is proportional to the difference in temperature from the outer wall to the surroundings, which can be modelled by
dθ
dr
r=2
= 0.083(θ(2) − 20)
Solve the equation with this boundary condition.
8. The governing equation for a projectile shot vertically upward is
m
d
2y
dt2
= −mg − C|V |V, y(0) = 0, y′
(0) = V0
where m is the mass of the projectile (kg), y(t) is the height (m), g is the acceleration of gravity (9.80665m/s2), C is an aerodynamic drag parameter, and V = dy/dt is 5 the velocity. For m = 10kg, C = 0.1Ns2/m2 , and V0 = 500.0m/s, calculate (a)
the maximum height attained by the projectile, (b) the time required to reach the maximum height, and (c) the time required to return to the original elevation.
9. A cable hung at its two ends as shown in Fig. 1 by its own weight will have a catenary shape described by the equation
y =
Tx
w
cosh wx
Tx
− 1
(1)
where w is the weight per unit length and Tx is the horizontal, x-component of the tension of the cable. Equation (1) is for the case when both w and Tx are constant throughout the cable. In fact, Equation (1) is the solution of the differential equation
d
2y
dx2
=
w
Tx
ds
dx =
w
Tx
"
1 +
dy
dx2
#1/2
(2)
where s is a variable along the length of the cable.
(a) Applying the Runge-Kutta method to solve the equation (2), where w = 0.12 kN/m, xA = yA = 0, xB = 200 m, and yB = 50 m. Let the initial conditions be y = 0 at x = xA = 0, iterate Tx value until yB is within 99.9% of 50 m.
(b) How could the problem be solved if xA = 100 m and yA = 25 m by application of the Runge-Kutta method (w remains equal to 0.12 kN/m)?
(c) Making a discussion about numerical results with different mesh steps and with other methods if possible.
10. The equation
y
′ = 1 + y
2
, y(0) = 0
has the analytical solution y(t) = tan t. The tangent function is infinite at t = π/2.
(a) Write your program by using both Euler and Runge-Kutta method to solve this initial value problem between t = 0 and t = 1.6, and compare the results of the program with the analytical solution
(b) What is the behaviour of the Runge-Kutta fourth order method when used between t = 0 and t = 1.6 6
(c) Compare the results from your codes and make comments on the behaviour change with step sizes.
11. Given the following non-linear boundary value problem
y
′′ − 18y
2 = 0, y(1) = 1
3
, y(2) = 1
12
(a) Use the shooting method to approximate solution
(b) Use finite difference to approximate solution
(c) Plot the approximate solutions together with the exact solution y(t) = 1 3t 2 and discuss your results with both methods
12. Here is a typical steady-state heat flow problem. Consider a thin steel plate to be a 10 × 20 (cm)2 rectangle. If one side of the 10 cm edge is held at 1000C and the other three edges are held at 00C, what are the steady-state temperature at interior points?
We can state the problem mathematically in this way if we assume that heat flows only in the x and y directions:
Find u(x, y) (temperature) such that
∂
2u
∂x2
+
∂
2u
∂y2
= 0 (3)
with boundary conditions
u(x, 0) = 0
u(x, 10) = 0
u(0, y) = 0
u(20, y) = 100
We replace the differential equation by a difference equation
1
h2
[ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j ] = 0 (4)
which relates the temperature at the point (xi
, yj ) to the temperature at four neighbouring points, each the distance h away from (xi , yj ). An approximation of Equation
(3) results when we select a set of such points (these are often called as nodes) and find the solution to the set of difference equations that result.
(a) If we choose h = 5 cm , find the temperature at interior points.
(b) Write a program to calculate the temperature distribution on interior points with h = 2.5, h = 0.25, h = 0.025 and h = 0.0025 cm. Discuss your solutions and examine the effect of grid size h.
(c) Modified the difference equation (4) so that it permits to solve the equation
∂
2u
∂x2
+
∂
2u
∂y2
= xy(x − 2)(y − 2)
on the region
0 ≤ x ≤ 2, 0 ≤ y ≤ 2
with boundary condition u = 0 on all boundaries except for y = 0, where u = 1.0.
Write and run the program with different grid sizes h and discuss your numerical results.
13. If a anti level beam of length L, which bends due to a uniform load of w lb/ft, is also subject to an axial force P at its free end (see Figure 2). The equation of its elastic curve is
EI d
2y
dx2
= P y −
wx2
2
(5)
For this equation, the origin O has been taken at the free end. I is the moment of inertia; here I =
bh3
12
. At the point x = L,
dy
dx = 0 and at x = 0, y = 0.
P
O
X
Y
L
(a) Solve the boundary value problem by any numerical method for a wooden beam, 2 in × 4 in × 10 ft, where E = 12 × 105 lb/in2
. Find y versus x when the beam has the 4 in dimension vertical with w = 25 lb/ft and a tension force of P = 500 lb.
(b) And also solve for the deflections if the beam is turned so that the 4 in dimension is horizontal.
(c) If the factor of the radius of curvature has been taken into account, Equation (5) becomes
EI
[1 + (y
′)
2]
3/2
d
2y
dx2
= P y −
wx2
2
(6)
If the defection of the beam is small, the difference is negligible, but in some cases this is not true. Furthermore, if there is considerable bending of the beam, the horizontal distance from the origin to the wall is less than L, the original length of the beam. Solve Equation (6) with the same conditions, and determine by how much the deflections differ from those previous calculations.
14. Consider a laterally insulated metal bar of length L = 1, which satisfies the heat equation
∂u
∂t =
∂
2u
∂x2
Suppose that the ends of the bar are kept at temperature u = 0oC and the temperature in the bar at t = 0, u(0, x) = sin πx.
(a) Using an explicit method with different meshes h to solve the equation and find the numerical solutions, then plot the solutions for t = 0, t = 1, t = 2, t = 3, t = 4 and t = 5 in the same figure.
(b) Using an implicit method with different meshes h to solve the equation and find the numerical solutions, then plot the solutions for t = 0, t = 1, t = 2, t = 3, t = 4 and t = 5 in the same figure.
(c) Using Crank-Nicolson method with different meshes h to solve the equation and find the numerical solutions, then plot the solutions for t = 0, t = 1, t = 2, t = 3, t = 4 and t = 5 in the same figure.
(d) Making your discussion for these three methods with different size of h The exact solution for this problem is
u(t, x) = sin(πx)e
−π
2
t
which can be used to compare with your numerical solutions.
15. Condon and Odishaw (1967) discuss Duffing equation for the flux φ in a transformer. This nonlinear differential equation is
d
2φ
dt2
+ ω
2
0φ + bφ3 =
ω
N
E cos(ωt), φ(0) = φ
′
(0) = 0
In this equation , E cos ωt is the sinusoidal source voltage and N is the number of turns in the primary winding, where ω0 and b are parameters of the transformer design. Make a plot of φ versus t (and compare to the source voltage) if E = 165, ω = 120π,
N = 600, ω
2
0 = 83 and b = 0.14. For approximate calculations, the nonlinear term bφ3 is sometimes neglected. Evaluate your result to determine whether this makes a significant error in the results. (Condon E. and Odishaw H., (1967) Handbook of Physics, New York, McGraw-Hill)
16. An 80kg paratrooper is dropped from an airplane at a height of 600m. After 5 second the chute opens. T