1
10.18 Write a user-defined MATLAB function that solves a first-order ODE by using the second-order
Adams–Bashforth method. The function should use the modified Euler method for calculating the solution
at the second point, and Eq. (10.95) for calculating the solution at the rest of the points. For function name
and arguments use [x,y]=odeAdams2(ODE,a,b,h,yINI). The input argument ODE is a name for
dy
the function that calculates ------ . It is a dummy name for the function that is imported into odeAdams2.
dx
The arguments a and b define the domain of the solution, h is the step size, and yINI is the initial value.
The output arguments, x and y, are vectors with the x and y coordinates of the solution.
Use the function odeAdams2 to solve the ODE in Problem 10.3. Write a MATLAB program in a
script file that solves the ODE by using h = 0.1 . The program should also plot the exact solution (given in
Problem 10.3) and the numerical solution (both in the same figure).
Solution
The listing of the user-defined function odeAdams2 is:
function [x, y] = odeAdams2(ODE,a,b,h,yINI)
% odeMIDPOINT solves a first order initial value ODE using the
% second-order Abam-Bashforth method.
% Input variables:
% ODE Name of a function file that calculates dy/dx.
% a The first value of x.
% b The last value of x.
% h Step size.
% yINI The value of the solution y at the first point (initial value.
% Output variable:
% x A vector with the x coordinate of the solution points.
% y A vector with the y coordinate of the solution points.
x(1) = a; y(1) = yINI;
N = (b-a)/h;
x(2) = x(1) + h;
K1=ODE(x(1),y(1));
K2=ODE(x(2),y(1)+K1*h);
y(2) = y(1) + (K1+K2)*h/2;
Excerpts from this work may be reproduced by instructors for distribution on a not-for-profit basis
for testing or instructional purposes only to students enrolled in courses for which the textbook
has been adopted. Any other reproduction or translation of this work beyond that permitted by
Sections 107 or 108 of the 1976 United States Copyright Act without the permission of the
copyright owner is unlawful.
2
for i = 2:N
x(i+1) = x(i) + h;
y(i+1) = y(i) + (3*ODE(x(i),y(i))-ODE(x(i-1),y(i-1)))*h/2;
end
The the user-defined function odeAdams2 is used in following program (script file) to solve Problem 8.3:
clear, clc
ODEHW10_3 = @ (t,y) y+t^3;
a=0; b=1.5; yINI=1;
h=0.1;
[x1, y1] = odeAdams2(ODEHW10_3,a,b,h,yINI);
fplot('7*exp(t)-t^3-3*t^2-6*t-6',[0 1.5])
hold on
plot(x1,y1,'*r')
hold off
xlabel('x'); ylabel('y')
legend('Exact','Numerical',2)
disp('Numerical solution with h=0.1')
Answer=[x1;y1]
When the program is executed the numerical solution is displayed in the Command Window and the fol-
lowing plot is displayed in the Figure Window:
Numerical solution with h=0.1
Answer =
Columns 1 through 7
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000
1.0000 1.1051 1.2210 1.3500 1.4951 1.6601 1.8499
Columns 8 through 14
0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000
2.0706 2.3293 2.6348 2.9973 3.4287 3.9428 4.5555
Columns 15 through 16
1.4000 1.5000
Excerpts from this work may be reproduced by instructors for distribution on a not-for-profit basis
for testing or instructional purposes only to students enrolled in courses for which the textbook
has been adopted. Any other reproduction or translation of this work beyond that permitted by
Sections 107 or 108 of the 1976 United States Copyright Act without the permission of the
copyright owner is unlawful.
3
5.2848 6.1515
7
Exact
Numerical
6
y
3
1
0 0.5 1 1.5
x
Excerpts from this work may be reproduced by instructors for distribution on a not-for-profit basis
for testing or instructional purposes only to students enrolled in courses for which the textbook
has been adopted. Any other reproduction or translation of this work beyond that permitted by
Sections 107 or 108 of the 1976 United States Copyright Act without the permission of the
copyright owner is unlawful.
