MATH3242.03/COSC3133.03

FW03

Solutions to Assignment 5

1. Let y2j,1 be the result of applying the modified midpoint method with h2j = H/2j.

Then we have

y(x+H) = y2j,1 + c1h2j2 + c2 h2j4+ c3 h2j6 + ... (1)

y(x+H) = y2j-2,1 + c1h2j-22 + c2h2j-24 + c3 h2j-26 + ... (2)

Thus evaluating (2j)2 (1) - (2j-2)2 (2) to eliminate the term in H2 we get

((2j)2 - (2j-2)2) y(x+H) = (2j)2 y2j,1 - (2j-2)2 y2j-2,1 + c2[H4/(2j)2 - H4/(2j-2)2] + O(H6)

so that defining

this becomes

y(x+H) = y2j,2 - c2 H4/(2j)2(2j-2)2 + O(H6) (3)

Similarly

y(x+H) = y2j-2,2 - c2 H4/(2j-2)2(2j-4)2 + O(H6) (4)

Now eliminating the H4 term between (3) and (4) yields

((2j)2 - (2j-4)2) y(x+H) = (2j)2 y2j,2 - (2j-4)2 y2j-2,2 +O(H6)

so that defining

we have

y(x+H) = y2j,3 + O(H6)

2. Below are the calculations involved in the modified midpoint method.

Here f(x,y) = 2x - 2w

 n=2, h=0.1 n=4, h=0.05 n=6, h=0.0333 x w x w x w 0. 1. 0. 1. 0. 1. 0.15 0.7 0.075 0.85 0.05 0.9 0.3 0.67 0.15 0.7675 0.1 0.83 y2,1 0.6295 0.225 0.66475 0.15 0.754 0.3 0.635575 0.2 0.7092 y4,1 0.624994 0.25 0.65216 0.3 0.628768 y6,1 0.6240256

Using the formula in from question 1 of the text we get the following extrapolation table

 2j y2j,1 y2j,2 y2j,3 2 0.6295 4 0.624994 0.623492 6 0.6240256 0.623251 0.623221
The exact result is 0.623217 so that y6,3 is accurate to 5 significant figures.

3. (a) y´ = -15y , 0 < x < 0.5, y(0) = 1. {exact solution is y = e-15x}

Euler method: wi+1 = wi - 15hwi = (1 - 15h)wi

Backward Euler method: wi+1 = wi - 15hwi+1 or wi+1 = wi/(1 + 15h)

Below are the results of the calculations with both h = 0.1 and 0.05

 x Euler Backward Euler Backward Exact 0. 1. 1. 1. 1. 1. 0.05 0.25 0.571429 0.472367 0.1 -0.5 0.4 0.0625 0.326531 0.223130 0.15 0.015625 0.186589 0.105399 0.2 0.25 0.16 0.003906 0.106622 0.049787 0.25 0.000977 0.060927 0.023518 0.3 -0.125 0.064 0.000244 0.034815 0.011109 0.35 0.000061 0.019895 0.005248 0.4 0.0625 0.0256 0.000015 0.011368 0.002479 0.45 0.000004 0.006496 0.001171 0.5 -0.03125 0.01024 0.000001 0.003712 0.000553

Note when h = 0.1 the Euler method alternates in sign. The Backward Euler method has the correct behaviour but is not very accurate for this value of h. Halving h means that both method have the right general behaviour but accuracy is still low.

(b) y´ = -15y + 225x, 0 < x < 0.5, y(0) = 0 {exact solution is y = e-15x + 15x - 1}

Euler method: wi+1 = wi + h(-15wi + 225xi) = (1 - 15h)wi + 225hxi

Backward Euler method: wi+1 = wi + h(-15wi+1 + 225hxi+1)

or wi+1 = (wi + 225hxi+1)/(1 + 15h)

Below are the calculations:

 x Euler Backward Exact 0. 0. 0. 0. 0.1 0. 0.9 0.7231 0.2 2.2500 2.1600 2.0479 0.3 3.3750 3.5640 3.5111 0.4 5.0625 5.0256 5.0025 0.5 6.4699 6.5102 6.5006

Here we see that both methods give good results for larger x since the exponential part decays quickly.

4. Solving x3 - 2z = 2 for x yields x = (2z + 2)1/3

Solving 5y2 - x3 = 7 for y yields y = ((x3 + 7)/5)1/2

Solving y2z = 1 for z yields z = y-2

These functions are continuous on the domain {1.3 < x < 1.6, 1.3 < y < 1.7, 0.3 < z < 0.8}. Furthermore on this domain (2z + 2)1/3 is a monotonically increasing function that lies in the interval (1.375, 1.533), ((x3 + 7)/5)1/2 is also a monotonically increasing function that lies in the interval (1.356, 1.490) while y-2 is a monotonically decreasing function that lies in the interval (0.346, 0.592). Thus the domain is mapped into itself by these functions.

The Jacobian for this system of equations has non-zero elements:

J13 = (2/3)(2z + 2)-2/3 < 0.35285 on the domain;

J21 = (3x2/10)((x3 + 7)/5)-1/2 = 0.3(5x4/(x3 + 7))1/2 < 0.51554;

J32 = -2y-3 < 0.91033.

Thus the infinity norm of J on the domain is max{0.35258, 0.51554, 0.91033} = 0.91033 < 1. Since the partial derivatives are continuous over the domain and the norm of the Jacobian is strictly less than unity there, the system has a unique fixed point in this domain.

The program and results for fixed-point iteration are given in a separate file.

5. Applying the linear shooting method to the second-order boundary-value problem

y´´ - 2xy´ + sin(x)y = (1-x)2, 0 < x < 2, y(0) = 1, y(2) = 1.5

we must solve the following two second-order initial-value problems for 0 ² x ² 2:

y1´´-2xy1´ + sin(x)y1 = (1-x)2, y1(0) = 1, y1´(0) = 0

and

y2´´- 2xy2´ + sin(x)y2 = 0, y2(0) = 0, y2 ´(0) = 1

Converting these into a system of first-order initial-value problems the first equation becomes:

u1´ = u2; u1(0) = 1

u2´ = 2xu2 - sin(x)u1 + (1-x)2; u2(0) = 0

and the second:

u3´ = u4; u3(0) = 0

u4´ = 2xu4 - sin(x)u3; u4(0) = 1

The required solution is

y(x) = y1(x) + (1.5 - y1(2)) y2(x)/y2(2)

The program and results are given in a separate file.

6. This is a nonlinear problem so we solve the initial value problem

y´ ´ = 2y y´ - xy2 + x cos(x), y(-1) = 1, y´ (-1) = t

We must also construct the differential equation for the z, the derivative of y with respect to the initial slope t. Since

f(x,y,y´) = 2y y´ - xy2 + x cos(x)

the partial derivative of f with respect to y is 2y´ - 2xy

and the partial derivative of f with respect to y´ is 2y. Thus the initial vaule problem for z is

z´ ´ = (2y´ - 2xy)z + 2yz´, z(-1) = 0, z´ (-1) = 1

We can convert these into a system of first-order equations with u1 = y and u3 = z

u1´ = u2, u1(-1) = 1

u2´ = 2u1u2 - xu12 + xcos(x), u2(-1) = t

u3´ = u4, u3(-1) = 0

u4´ = 2(u1u4 - xu1u3 + u2u3), u4(-1) = 1

The recursion relation for the slope t is tk+1 = tk - (u1(1,tk) + 1)/u3(1,tk)= tk - (y(1,tk) + 1)/z(1,tk)

7. In the notation of the notes p(x) = 2x , q(x) = -sin(x) and r(x) = (1-x)2. The program and results are given in a separate file.

5