Solving Adams-Bashforth Integration for ODE

This document explains the Adams-Bashforth integration method for solving ordinary differential equations (ODEs). It provides the derivation of the method using polynomial interpretation and the interpolation formula. The document also discusses the iteration equation and provides examples. The subject is MEC3456 and the document type is LAB 07.

Running head: MEC3456 LAB 07
MEC3456 LAB 07
Name of the Student
Name of the University
Author Note

2MEC3456 LAB 07
Question 1:
The Adams-Bashforth integration for finding the solution of the ODE dy/dt = f(t,y) are given
by the following integral update formula
yn+ 1= yn +∫
ti
ti+1
f ( t , y ) dt (1)
Where, yn+ 1 = y(tn+1 ¿
The method can be derived using polynomial interpretation as given below.
From equation (1) using the fundamental theorem of calculus we get
y(tn+1 ¿ = y(tn ¿ + ∫
tn
tn+1
y' ( t ) dt
Here, let A = ∫
tn
t n+1
y' ( t ) dt=∫
tn
t n+ 1
f (t , y ( t ) )dt
Now, for calculating the value of A the value of the interpolating polynomial P(t) as an
approximation f(t,y(t)) is used. Hence, the value of the interpolating polynomial in the
Lagrange form is
L(x) := ∑
j=0
k
yjlj(x ), where the expression of lj ( x )= ∏0 ≤m ≤ k∧m ≠ j
x −xm
xj−xm
Hence, the interpolation formula will be
P(t) = f ( tn , yn )∗t−tn−1
tn−tn−1
+f ( tn−1 , yn−1 )∗t−c
tn−1−tn
Hence, A = ∫
tn
tn+1
f ( t , y ( t ) ) dt ∫
tn
tn+1 f ( tn , yn ) ( t−tn−1 )
tn −tn−1
+ f ( tn−1 , yn−1 )∗t−tn
tn−1−tn
dt
Now, after integration and simplification gives,

3MEC3456 LAB 07
(½)(tn+tn +1 ¿ ( f ( tn , yn ) −f ( tn−1 , yn−1 ) )−tn−1 f ( tn , yn ) +tn f (tn−1 , yn−1)
= (½)( tn+ tn +1−2 tn−1 ¿ f ( tn , yn) + (½)(2tn−tn+1−t n ¿ f (tn−1 , yn−1 )
Now, tn−1 , c and tn+1 are equally spaced hence it can be considered that
tn−¿ tn−1 = tn+1−¿ tn = h. Hence, the value of A will be
A = (3/2)h*f( tn , yn ¿ – (½)h*f(tn−1 , yn−1 ¿
Hence, by Adams-Bashforth method the iteration equation will be this
y(tn+1 ¿ = y(tn ¿+ ( 3 h
2 )f ( tn , y n ) −( h
2 )f ( tn−1 , yn−1 )
Question 2:
2a:
The given second order differential equation is
( d
dx ) ( x3∗du
dx ) + λxu=0 (1)
The second order equation is defined inside the interval [1, 2] having boundary conditions
u(1) = 0 and u(2) = 0.
Hence, the above equation can be rewritten as
3 x2∗( du
dx )+ x3∗d2 u
d x2 + λx u= 0 (2)
Now, approximating (2) by finite difference method at an arbitrary point xi
3 xi
2∗ui+1 – ui
h + xi
3∗ui+1 – 2∗ui+ui−1
h2 + λ xi ui =0
Where, ui = u(xi) and holds true for all i.

4MEC3456 LAB 07
Q2b:
Now, by relaxation method the equation is first converted to two unknowns of equation
du
dx = y
x^3* d2 u
d x2 = - λxy− λu −3 x2∗y
dy/dx ¿ ¿- y− λu −3 x2∗y)/ x^3
In this method at first discretization by the method of finite differences is performed with the
mesh spacing of h= 0.25 for approximating the eigenvalue problem. The given boundary
conditions are u(1) = 0, u(2) = 0. An initial guess of λ = 0 is assumed and the solution is
improved iteratively. The first order differential equations are approximated by trapezoidal
rule here.
Q2c:
Now, the second order differential equation is converted to two first order equations.
du
dx = y
x^3* d2 u
d x2 = - λxy− λu −3 x2∗y
d2 u
d x2 =¿- λxy− λu −3 x2∗y)/ x^3
Now, u(2) = 0 and u(1) = 0 can be initialized to solve the system numerically by evaluating
the next iteration u(3). Here, the step size need to be h= 1.

5MEC3456 LAB 07
Here, any boundary condition can be chosen of u(1) and u’(1) as if u(x) is a solution of the
differential equation, then Cu(x) (C is any constant) is also a solution of the differential
equation. Now, for convenience we have chosen u(1) = 0 and u’(1) = 1.
Question 3a:
The uRHS value for x=2 for the differential equation is calculated by the following
MATLAB code.
MATLAB code:
function uRHS = RHS_Value(lambda)
% Calculates the RHS value of u given a guess for the eigenvalue lambda.
xspan = [1;2];
intcon = [0;1];
[x,u] = ode45(@odefunc,xspan,intcon);
uval = u(:,1);
% plot(x,uval,x,u(:,2))
% xlabel('x range in [1 2]')
% ylabel('solution u(x) and derivative du/dx')
% legend('u(x)','du/dx','Location','best')
uRHS = uval(end);

6MEC3456 LAB 07
function rhs = odefunc(x,u)
rhs = [u(2);(-lambda.*x.*u(2) - lambda.*u(1) - 3.*x.^2.*u(2))./(x.^3)];
end
end
Output:
uRHS =
0.2869
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x range in [1 2]
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
solution u(x) and derivative du/dx
u(x)
du/dx
Question 3b:

7MEC3456 LAB 07
Now, for λ = [0,1000] the values of u(2) are plotted by the following code.
MATLAB code:
lambdarange = 0:1:1000;
for i=1:length(lambdarange)
uRHS(i) = RHS_Value(lambdarange(i));
end
figure(1)
plot(lambdarange,uRHS)
xlabel('\lambda')
ylabel('u(2)')
title('u value at x=2 for \lambda in [0,1000]')
Plot:

8MEC3456 LAB 07
0 100 200 300 400 500 600 700 800 900 1000
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
u(2)
u value at x=2 for in [0,1000]
Question 3c:
Given, differential equation
x^3* d2 u
d x2 = - λx∗du
dx − λu −3 x2∗du
dx
The above equation is of the form of Euler differential equation and the roots of the equation
that is required to be solved is
r(r-1)(r-2) + (3r*(r-1) + λr) + λ = 0
The solution has the form
u ( x ) =( 1
x )∗sin ( nπ∗ln ( x )
ln ( 2 ) )
Hence, solving for r gives,

9MEC3456 LAB 07
r1 = -1, r2 = ½− ( 1−4∗λ )
1
2
2
, r3 = ( 1 – 4∗λ )
1
2
2 +½
Now, for real r the eigenvalues should be
λn=1+ ( nπ
ln 2 ) 1
2
MATLAB code:
% Q3c - find eigenvalues
n = 1:5;
lambda = 1 + ((n.*pi)./(log(2))).^(1/2);
sprintf('The first five eigenvalues are \n %.5f \n %.5f \n %.5f \n %.5f \n
%.5f',lambda(1),lambda(2),lambda(3),lambda(4),lambda(5))
Output:
'The first five eigenvalues are
3.12893
4.01077
4.68742
5.25787
5.76044'
Question 3d:

10MEC3456 LAB 07
Now, for first and fourth eigenvalues the analytical solution is plotted by the following
MATLAB code.
MATLAB code:
lamdafirst = lambda(1); lamdafourth = lambda(4);
x = 1:0.01:2;
uvalfirst = x.^(-1).*sin(lamdafirst*pi*log(x)./log(2));
uvalfourth = x.^(-1).*sin(lamdafourth*pi*log(x)./log(2));
plot(x,uvalfirst)
hold on
plot(x,uvalfourth)
ylim([-1 1])
xlabel('x value')
ylabel('u(x) for two eigenvalues')
legend('Analytical u(x) for first eigenvalue','anlaytical u(x) for second eigenvalue')
Plot:

11MEC3456 LAB 07
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x value
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
u(x) for two eigenvalues
Analytical u(x) for first eigenvalue
anlaytical u(x) for second eigenvalue
Question 4:
Now, the differential equation is solved using relaxation method in MATLAB. The first five
eigenvalues λ1, λ2, λ3, λ4 and λ5 are listed and the corresponding eigenvectors are obtained.
Furthermore the eigenvectors for λ1 and λ4 are plotted.
The given differential equation is
3 x2∗( du
dx )+ x3∗d2 u
d x2 + λx u= 0
Now, dividing both sides by x gives,
x^2* d2 u
d x2 + 3x*( du
dx )+ λ u=0
Employing finite difference approach

12MEC3456 LAB 07
3 xi∗ui+1 – ui
h + xi
2∗ui+1 – 2∗ui+ui−1
h2 + λ ui=0
Now, the boundary conditions are u(1) = u(2) = 0
Hence, by relaxation approach with n=20 unknowns y1 = 0 and y22 = 0.
Where, y2 to y21 are 20 unknown u values in the range x=[1,2] and y1 = u(1) = 0 and y22 =
u(2) = 0.
Hence, the general equations in terms of internal nodes will be
3 x2∗ y3 – y2
h + x2
2∗y3 – 2∗y2 + y1
h2 + λ y2=0
3 x3∗y4 – y3
h + x3
2∗y 4 – 2∗y3 + y2
h2 + λ y3=0
….
The last equation is
3 x2 1∗y22 – y21
h + x2 1
2 ∗ y22 – 2∗y2 1+ y20
h2 + λ y2 1=0
Now, putting boundary conditions and separating λ term in right side
3 x2∗ y3 – y2
h + x2
2∗y3 – 2∗y2
h2 =−λ y2
3 x3∗y4 – y3
h + x3
2∗y 4 – 2∗y3 + y2
h2 =−λ y3
….

13MEC3456 LAB 07
3 x21∗– y21
h + x21
2 ∗– 2∗y21 + y20
h2 =−λ y21
Hence, in matrix form
[ −3 x 2
h − 2
h2 x 22 3 x 2
h + x 22
h2 … ..0
x 32
h2
−3 x 3
h − 2
h2 x 32 3 x 3
h + x 32
h2
0 … 0 x21
2
h2
−3 x 21
h − 2 x 2 12
h2
] [ y 2
…
y 21 ]=−λ [ y 2
…
y 21 ]
MATLAB code:
n=20;
x = linspace(1,2,n+2);
h = x(2) - x(1);
A = zeros(n);
for i=2:n-1
j=i-1;
A(i,j) = (x(i+1)^2)/h^2;
A(i,j+1) = -3*x(i+1)/h - 2*(x(i+1)^2)/h^2;
A(i,j+2) = (3*x(i+1))/h + (x(i+1)^2)/h^2;
end

14MEC3456 LAB 07
% defining 1st and last row of matrix A
A(1,1) = -3*x(2)/h - (2*x(2)^2)/h^2; A(1,2) = 3*x(2)/h + x(2)^2/h^2;
A(20,19) = (x(21)^2)/h^2; A(20,20)= - 3*x(21)/h - 2*x(21)^2/h^2;
% Representing A in Ay = lambda*y form
A = -A;
[V,lambda]= eig(A);
for i=1:length(A)
eigenvals(i,1) = lambda(i,i);
end
eigenvalseq = sort(eigenvals);
fprintf('The first five eigen values are \n%.4f \n%.4f \n%.4f \n%.4f \n%.4f \
n',eigenvalseq(1),eigenvalseq(2),eigenvalseq(3),eigenvalseq(4),eigenvalseq(5))
for i=1:length(A)
if eigenvals(i,1) == eigenvalseq(1) % extracting position of first eigenvalue in sorted
manner
pos1 = i;
elseif eigenvals(i,1) == eigenvalseq(4) % extracting position of fourth eigenvalue in sorted
manner
pos4 = i;

15MEC3456 LAB 07
else
end
end
solfirst = V(:,pos1);
solfourth = V(:,pos4);
plot(x(2:n+1)',solfirst,x(2:n+1)',solfourth)
grid on
title('Eigen solution for first and fourth eigenvalue')
legend('Solution vector for first eigenvalue','Solution vector for fourth
eigenvalue','Location','best')
xlabel('x value in [1,2]')
ylabel('u(x) for two different eigenvalues')
Output for n=20:
The first five eigen values are
22.5105
86.6105
191.6550
334.9829
512.9095

16MEC3456 LAB 07
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x value in [1,2]
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
u(x) for two different eigenvalues
Eigen solution for first and fourth eigenvalue
Solution vector for first eigenvalue
Solution vector for fourth eigenvalue
Now, for n=100 the value of n in the first line of script is changed.
Output for n=100:
The first five eigen values are
-0.0000
22.0532
85.2128
190.4011
337.5015

17MEC3456 LAB 07
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x value in [1,2]
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
u(x) for two different eigenvalues
Eigen solution for first and fourth eigenvalue
Solution vector for first eigenvalue
Solution vector for fourth eigenvalue