FW03
Solutions to Assignment 5
1.
Let y_{2j,1} be the result of applying the modified midpoint
method
with h_{2j} = H/2j.
Then we have
y(x+H) = y_{2j,1} + c_{1}h_{2j}^{2} + c_{2 }h_{2j}^{4}+ c_{3 }h_{2j}^{6 }+ ... (1)
y(x+H) = y_{2j2,1} + c_{1}h_{2j2}^{2 }+ c_{2}h_{2j2}^{4} + c_{3 }h_{2j2}^{6} + ... (2)
Thus evaluating (2j)^{2} (1)  (2j2)^{2} (2) to eliminate the term in H^{2} we^{ }get
((2j)^{2}  (2j2)^{2}) y(x+H) = (2j)^{2} y_{2j,1}  (2j2)^{2} y_{2j2,1} + c_{2}[H^{4}/(2j)^{2 } H^{4}/(2j2)^{2}] + O(H^{6})
so that defining
this becomes
y(x+H) = y_{2j,2}  c_{2 }H^{4}/(2j)^{2}(2j2)^{2} + O(H^{6}) (3)
Similarly
y(x+H) = y_{2j2,2}
 c_{2 }H^{4}/(2j2)^{2}(2j4)^{2
}+ O(H^{6})
(4)
Now eliminating the H^{4} term between (3) and (4) yields
((2j)^{2}  (2j4)^{2}) y(x+H) = (2j)^{2} y_{2j,2}  (2j4)^{2} y_{2j2,2} +O(H^{6})
so that defining
we have
y(x+H) = y_{2j,3} + O(H^{6})
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 
y_{2,1}  0.6295  0.225  0.66475  0.15  0.754 
0.3  0.635575  0.2  0.7092  
y_{4,1}  0.624994  0.25  0.65216  
0.3  0.628768  
y_{6,1}  0.6240256 
Using the formula in from question 1 of the text we get the following extrapolation table
2j  y_{2j,1}  y_{2j,2}  y_{2j,3} 
2  0.6295  
4  0.624994  0.623492  
6  0.6240256  0.623251  0.623221 
3. (a) y´ = 15y , 0 < x < 0.5, y(0) = 1. {exact solution is y = e^{15x}}
Euler method: w_{i+1} = w_{i}  15hw_{i} = (1  15h)w_{i}
Backward Euler method: w_{i+1} = w_{i}  15hw_{i+1} or w_{i+1} = w_{i}/(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: w_{i+1} = w_{i} + h(15w_{i} + 225x_{i}) = (1  15h)w_{i} + 225hx_{i}
Backward Euler method: w_{i+1} = w_{i} + h(15w_{i+1} + 225hx_{i+1})
or w_{i+1} = (w_{i} + 225hx_{i+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 x^{3}  2z = 2 for x yields x = (2z + 2)^{1/3}
Solving 5y^{2}  x^{3} = 7 for y yields y = ((x^{3} + 7)/5)^{1/2}
Solving y^{2}z = 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), ((x^{3} + 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 nonzero elements:
J_{13} = (2/3)(2z + 2)^{2/3} < 0.35285 on the domain;
J_{21} = (3x^{2}/10)((x^{3} + 7)/5)^{1/2} = 0.3(5x^{4}/(x^{3} + 7))^{1/2} < 0.51554;
J_{32} = 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 fixedpoint iteration are given in a separate file.
5. Applying the linear shooting method to the secondorder boundaryvalue problem
y´´  2xy´ + sin(x)y = (1x)^{2}, 0 < x < 2, y(0) = 1, y(2) = 1.5
we must solve the following two secondorder initialvalue problems for 0 ² x ² 2:
y_{1}´´2xy_{1}´ + sin(x)y_{1} = (1x)^{2}, y_{1}(0) = 1, y_{1}´(0) = 0
and
y_{2}´´ 2xy_{2}´ + sin(x)y_{2} = 0, y_{2}(0) = 0, y_{2 }´(0) = 1
Converting these into a system of firstorder initialvalue problems the first equation becomes:
u_{1}´ = u_{2}; u_{1}(0) = 1
u_{2}´ = 2xu_{2}  sin(x)u_{1} + (1x)^{2}; u_{2}(0) = 0
and the second:
u_{3}´ = u_{4}; u_{3}(0) = 0
u_{4}´ = 2xu_{4}  sin(x)u_{3}; u_{4}(0) = 1
The required solution is
y(x) = y_{1}(x) + (1.5  y_{1}(2)) y_{2}(x)/y_{2}(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´  xy^{2} + 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´  xy^{2} + 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 firstorder equations with u_{1} = y and u_{3} = z
u_{1}´ = u_{2, }u_{1}(1) = 1
u_{2}´ = 2u_{1}u_{2}  xu_{1}^{2} + xcos(x), u_{2}(1) = t
u_{3}´ = u_{4}, u_{3}(1) = 0
u_{4}´ = 2(u_{1}u_{4}  xu_{1}u_{3} + u_{2}u_{3}), u_{4}(1) = 1
The
recursion relation for the slope t is t_{k+1 }= t_{k}

(u_{1}(1,t_{k}) + 1)/u_{3}(1,t_{k})= t_{k}

(y(1,t_{k}) + 1)/z(1,t_{k})
7.
In the notation of the notes p(x) = 2x , q(x) = sin(x) and r(x) = (1x)^{2}.
The program and results are given in a separate file.