Chapter 4 Exercises

4.7. Chapter 4 Exercises#

You should try the following exercise questions first, then check with the answers.

Exercise 4.1

Use a G-N interpolation formula to derive the following multistep formulae:

$ y_{j+2} = y_j + \frac{h}{3}(f_{j+2} + 4f_{j+1} + f_j) $ [ G-N forward 0 to 2 ]
$ y_{j+2} = y_{j+1} + \frac{h}{12}(5f_{j+2} + 8f_{j+1} - f_j) $ [ G-N forward 1 to 2 ]
$ y_{j+3} = y_{j-1} + \frac{h}{3}(8f_{j+2} - 4f_{j+1} + 8f_j) $ [ G-N forward -1 to 3 ]
$ y_j = y_{j-1} + \frac{h}{12}(5f_j + 8f_{j-1} - f_{j-2}) $ [ G-N backward 0 to -1 ]

Solution to

handwritten full solution for Exercise 3.1

Exercise 4.2

The family of Multistep explicit Adams methods can be expanded as

\[\begin{aligned} y_{j+1} ~&=~ y_j + h(1 + \frac{1}{2}\nabla + \frac{5}{12}\nabla^2 + \frac{3}{8}\nabla^3 + \frac{251}{720}\nabla^4 + \dots)f_j \end{aligned}\]

Use this formula and the following difference table to solve the initial value problem:

\[\begin{aligned} y' ~&=~ -y - 2x, & y(0) ~&=~ -1, & \text{for $x = 0.1$ to $x = 0.5$ with $h = 0.1$.} \end{aligned}\]

Work to four decimal places. Compare your results with the analytical solution: $\,y(x) = -3e^{-x} - 2x + 2\,$. Comment on the accuracy of your results.

$j$	$x$	$y$	$f$	$\nabla f$	$\nabla^2 f$
0	0.0	-1.00000	1.0000
				…………
1	0.1	…………	…………		…………
				…………
2	0.2	…………	…………

Solution to

$y_1 = -0.9000$, $y_2 = -0.8450$, $y_3 = -0.8114$, $y_4 = -0.8019$, $y_5 = -0.8102$.

Exercise 4.3

The family of Multistep explicit Adams methods can be expanded as

(4.6)#\[\begin{aligned} y_{j+1} ~&=~ y_j + h(1 + \frac{1}{2}\nabla + \frac{5}{12}\nabla^2 + \frac{3}{8}\nabla^3 + \frac{251}{720}\nabla^4 + \dots)f_j \end{aligned}\]

Similarly, the family of Multistep implicit Adams methods can be expanded as

(4.7)#\[\begin{aligned} y_{j+1} ~&=~ y_j + h(1 - \frac{1}{2}\nabla - \frac{1}{12}\nabla^2 - \frac{1}{24}\nabla^3 - \frac{19}{720}\nabla^4 - \dots)f_{j+1} \label{eq:ch02:2} \end{aligned}\]

Use formula (4.6) as a predictor, formula (4.7) as a corrector and the following difference table to solve the initial value problem: $y' = -y -2x$, $y(0) = -1$, to advance the solution to $x = 0.4$. State results to four decimal places. Compare your result at each step with the corresponding solution obtained in question 2, above, and comment.\

$j$	$x$	$y$	$f$	$\nabla f$	$\nabla^2 f$	$\nabla^3 f$	$\nabla^4 f$	$\nabla^5 f$
0	0.0	-1.00000	1.0000
				…………
1	0.1	…………	…………		…………
				…………		…………
2	0.2	…………	…………		…………		…………
				…………		…………		…………
3	0.3	…………	…………		…………		…………
				…………		…………
4	0.4	…………	…………		…………
				…………
5	0.5	…………	…………

Solution to

$x = 0.1$, $y^p = -0.9000$, $y^c = -0.9150$, $y_{ex} = -0.9145$
$x = 0.2$, $y^p = -0.8577$, $y^c = -0.8566$, $y_{ex} = -0.8562$
$x = 0.3$, $y^p = -0.8227$, $y^c = -0.8228$, $y_{ex} = -0.8225$
$x = 0.4$, $y^p = -0.8112$, $y^c = -0.8113$, $y_{ex} = -0.8110$

excel solution
A video explanation

Exercise 4.4

Consider the initial–value problem

\[\begin{aligned} y' ~&=~ x^2 - y, & y(0) ~&=~ 1. \end{aligned}\]

Use the step-by-step method applied in Example 4.2, to find approximate solutions at $x = 0.4$ and $x = 0.5$ using the $4^\text{th}$ order ABM predictor-corrector method (listed on page 6 Chapter 2), with $\,h = 0.1\,$. Use the exact solution, $y_{ex} = 2 - 2x + x^2 - e^{-x}$, to calculate the starting values. Work to 7 decimal places. Note: you can use an Excel spreadsheet to carry out your calculations.
Modify Program 1 (Matlab and Python) to find approximate solutions in the range $x = 0$ to $x=0.5$, with $\,h = 0.1\,$. Use the exact solution, $y_{ex} = 2 - 2x + x^2 - e^{-x}$, to calculate the starting values.
Modify Program 2 (Matlab and Python) to check your results at $x = 0.4$ and $x = 0.5$. Run your program for smaller and larger step size values and comment on your results.
Modify Program 3 (Matlab only) - with the ode113 solver - to find solution and comment on your results. Try different tolerance values for AbsTol and RelTol and comment on your findings.

Solution to

Answer

$h = 0.1$,
$x = 0.4$, $y^p = 0.6896771$, $y^c = 0.6896803$, $y_{ex} = 0.6896800$
$x = 0.5$, $y^p = 0.6434670$, $y^c = 0.6434699$, $y_{ex} = 0.6434693$ $y_{ex} = 0.6434693$

excel solution

Answer

Case 1: RelTol=AbsTol=1.e-4

$x = 0.4$, $y = 0.6896797$, $y_{ex} = 0.6896800$, $|y_{ex} - y| = 0.00000023$
$x = 0.5$, $y = 0.6434691$, $y_{ex} = 0.6434693$, $|y_{ex} - y| = 0.00000021$

Matlab code

%Exercise 3.4 Question 4
%Scriptfile to solve a single 1st order ordinary differential equation, using
%a non-stiff ODE solver, called ode113, from the Matlab ODE solver routines
%based on Adam-Bashforth-Moulton methods.
function ch2_prog3
y0=1.0;  %set the initial value of the dependent variable y
%tspan=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]; %set specific output point
tspan=0:0.1:1;
options = odeset('RelTol',1e-4,'AbsTol',1e-4); %set tolerances
[t,y]=ode113(@f,tspan,y0,options); %call the ODE solver routine ode113


disp('  t        y          f(t)      y_ex     Abs_error ');


for i=1:length(t)  %in this case length(t)=11, from tspan
    yex(i)=yexact(t(i));  %calculating the exact solutions at each step
    abserror(i)=abs(yex(i)-y(i));  %calculating the absolute error
    fprintf('%5.2f %10.7f %10.7f %10.7f %11.8f\n', t(i),y(i),f(t(i),y(i)),yex(i),abserror(i));
end

plot(t,y(:,1),'LineWidth',2);
title(['Solution of dy/dt = -y+t+1, y0 = ' num2str(y0),', using ode113']);
xlabel('t');
ylabel('solution y');
%--------------------------------------------------------------------------
%Defining the ODE to be solved
function f=f(t,y)
f=t.^2-y;

% The exact solution
function yexact=yexact(x)
    yexact = 2-2*x+x^2-exp(-x);        

Case 2: RelTol=AbsTol=1.e-6

$x = 0.4$, $y = 0.6896799$, $y_{ex} = 0.6896800$, $|y_{ex} - y| = 0.00000001$
$x = 0.5$, $y = 0.6434693$, $y_{ex} = 0.6434693$, $|y_{ex} - y| = 0.00000001$

Matlab code

%Exercise 3.4 Question 4
%Scriptfile to solve a single 1st order ordinary differential equation, using
%a non-stiff ODE solver, called ode113, from the Matlab ODE solver routines
%based on Adam-Bashforth-Moulton methods.
function ch2_prog3
y0=1.0;  %set the initial value of the dependent variable y
%tspan=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]; %set specific output point
tspan=0:0.1:1;
options = odeset('RelTol',1e-6,'AbsTol',1e-6); %set tolerances
[t,y]=ode113(@f,tspan,y0,options); %call the ODE solver routine ode113


disp('  t        y          f(t)      y_ex     Abs_error ');


for i=1:length(t)  %in this case length(t)=11, from tspan
    yex(i)=yexact(t(i));  %calculating the exact solutions at each step
    abserror(i)=abs(yex(i)-y(i));  %calculating the absolute error
    fprintf('%5.2f %10.7f %10.7f %10.7f %11.8f\n', t(i),y(i),f(t(i),y(i)),yex(i),abserror(i));
end

plot(t,y(:,1),'LineWidth',2);
title(['Solution of dy/dt = -y+t+1, y0 = ' num2str(y0),', using ode113']);
xlabel('t');
ylabel('solution y');
%--------------------------------------------------------------------------
%Defining the ODE to be solved
function f=f(t,y)
f=t.^2-y;

% The exact solution
function yexact=yexact(x)
    yexact = 2-2*x+x^2-exp(-x);        