BRUNEL UNIVERSITY
ASSIGNMENT MATHS AND COMPUTATIONAL METHODS
2012
ME2602
Omar user
3/28/2012
ASSIGNMENT MATHS AND COMPUTATIONAL METHODS 2012
1. Newton Interpolation
a)
n=length(x);
%First two columns of divided difference table which are the given x and y
%values
DD(1:length(x),1) = x;
DD(1:length(y),2) = y;
%Below creates a table of coefficients b0, b1 etc. using divided
difference.
for j = 1:n-1;
for k = 1:n-j;
DD(k,j+2) = (DD(k+1,j+1)-DD(k,j+1))/(DD(k+j,1)-DD(k,1));
end
end
DD(6,1)=0;
DD(7,1)=0;
%This puts coefficients for each order in an array so they can be used by
%f(k)
Result for divided difference table from above code:
b0
b1
b2
b3
b4
=
=
=
=
=
y;
[DD(1,3),DD(2,3),DD(3,3),DD(4,3),DD(5,3)];
[DD(1,4),DD(2,4),DD(3,4),DD(4,4),DD(5,4)];
[DD(1,5),DD(2,5),DD(3,5),DD(4,5),DD(5,5)];
[DD(1,6),DD(2,6),DD(3,6),DD(4,6),DD(5,6)];
%Newton interpolation formula to calculate interpolations of 1st to 5th
orders.
%x(6) and x(7) must be called 0 in order for iterations to continue, in the
for loop.
x(6) = 0;
x(7) = 0;
%This is the application of the newton interpolation formula where x1 is
the point at which interpolation is occurring.
for k = 1:4;
f(k)=b0(k)+b1(k)*(x1-x(k))+b2(k)*(x1-x(k))*(x1-x(k+1))+b3(k)*(x1-x(k))*(x1x(k+1))*(x1-x(k+2))+b4(k)*(x1-x(k))*(x1-x(k+1))*(x1-x(k+2))*(x1-x(k+3));
end
%Prints requested results.
fprintf('Interpolating at x = 1/3\n')
fprintf('\n');
fprintf('First order interpolation gives: %g\n',f(1,4))
fprintf('Second order interpolation gives: %g\n',f(1,3))
fprintf('Third order interpolation gives: %g\n',f(1,2))
fprintf('Fourth order interpolation gives: %g\n',f(1,1))
b) Inserting the arrays below into the top of the previous script yields the following results on
Matlab.
Plotting the given function values
indicate the function could be even
due to the symmetry in the y axis.
Inserting the arrays below into the top of the previous script yields the following results on
Matlab.
% x and y values put into an array for each variable.
x = [-0.5,-1/6,0,1/6,0.5];
y = [2,1,0,1,2];
x1 = (1/3);
4th order newton extrapolation at x=1:
Matlab code for fourth order extrapolation (positioned below the previous function for the first
second third fourth and fifth order newton extrapolations):
%Newton extrapolation x=1.
xe=1;
for k = 1;
Ne(k)=b0(k)+b1(k)*(xe-x(k))+b2(k)*(xe-x(k))*(xe-x(k+1))+b3(k)*(xex(k))*(xe-x(k+1))*(xe-x(k+2))+b4(k)*(xe-x(k))*(xe-x(k+1))*(xe-x(k+2))*(xex(k+3));
end
fprintf('\n');
fprintf('Extrapolating at x = 1\n')
fprintf('Fourth order extrapolation gives: %g\n',Ne(1,1))
Result:
The extrapolation is not accurate as opposed to interpolating, as the solution lies outside the interval
used in the interpolation. This can cause large errors and give inaccurate solutions.
2.
a) Bisection method:
( )
For epsilon = 0.1
2(Xa)+1-tan(0.5*Xa)
Initial values where X(b)>X(a)
X(a)
-1.5708
-1.5708
-0.7854
-0.7854
-0.7854
-0.68722
X(b)
1.5708
0
0
-0.393
-0.589
-0.589
f(Xa)
-1.141593
-1.141593
-0.156583
-0.156583
-0.156583
-0.016641
f(Xb)
3.141593
1
1
0.413514
0.125249
0.125249
(Xa+Xb)/2
X(m)
0
-0.7854
-0.3927
-0.58905
-0.68722
-0.63814
f(Xm)
1
-0.1565828
0.4135142
0.1252494
-0.0166411
0.0540834
f(Xa) x f(Xm) f(Xa) *f(Xb)>0
-1.1415927
FALSE
0.1787537
TRUE
-0.0647492
FALSE
-0.0196119
FALSE
0.0026057
TRUE
-0.0009
FALSE
b) Newton Raphson:
For epsilon = 0.1
Initial approximation
2(X(0))+1-tan(0.5*X(0))
X(0)
f(X(0)
f'(X(0))
-1.5708 -1.142
1
-0.4292 0.3596 1.4762471
-0.67276 0.0041 1.4388656
2-0.5*(1/sec^2(0.5*X(0)))
X imp
-0.429
-0.673
-0.676
Ximp=X(0)-(f(X(0))/f(X(0))
c) Successive approximation method:
Ximp=g(X0)=0.5*tan(0.5*X(0))-0.5
For epsilon = 0.001
Initial approximation
X(0)
-1.571
-1.000
-0.773
-0.704
-0.684
-0.678
X imp
-1.000
-0.773
-0.704
-0.684
-0.678
-0.676
X imp -X(0)
0.5707963
0.2268488
0.0696225
0.0200141
0.005657
0.0015916
dg/dx
0.25
0.25
0.25
0.25
0.25
0.25
( ( ) ( )) ( ( ))
( ( )) ( ( ))
X(b)-X(a) ROOT x
3.141593
-0.733
1.570796
-0.733
0.785398
-0.679
0.392699
-0.678
0.19635
-0.676
0.098175
-0.676
3.
a)
( )
( )
Applying initial condition:
( )
( )
( )
b) Linear relationship:
Interpolating for lambda:
t(h)
0
1
2
3
4
5
6
7
8
9
t(h)
45
((
N(t)
10
18
30
42
70
120
205
333
550
880
N(t)
2258
( ) ( ( ))
((
Therefore:
ln(N(t))
2.30
2.89
3.40
3.74
4.25
4.79
5.32
5.81
6.31
6.78
ln(N(t))
45.59
( )
( ) )
t(h)^2
0
1
4
9
16
25
36
49
64
81
t(h)^2
285
( )) )
( ( ))))
t(h)ln(N(t))
0.00
2.89
6.80
11.21
16.99
23.94
31.94
40.66
50.48
61.02
t(h)ln(N(t))
245.93
t(h)
0
1
2
3
4
5
6
7
8
9
N(t)
10.00
16.39
26.88
44.06
72.23
118.41
194.12
318.24
521.72
855.30
N values calculated using
solution to differential
equation.
Plot of results
c)
The scenario makes sense as N (number of bacteria) stays almost constant or changes very little
after certain period of time. The graph below supports this where the constant A has been taken to
be 1.
Applying the condition that
(
Expressing
as partial fractions:
When
( )
:
(
When
( )
:
(
( )
Integrating:
)
( )
( )
( )
( )
a)
Rearranging for
( )
(
(
4.
)
( )
( )
(
)
:
( )
Expressing h in terms of M:
( )
(
)
)
(
b) Matlab code for iterating expression:
%Inputting variables for iterating.
%dt is the time step.
dt=0.14;
alpha = (2.3*10^-5);
xm = 0.1;
%M here is the number of points.
M= 20;
T0 = 300;
T1 = 200;
%Setting up a table to hold the values of the iterations
T(1,1:(M-1))=T0;
T(1:M,M)=T1;
T(2:M,1)=T0;
%Application of the iterative formula from part a to fill out the
constructed table above.
for j = 2:M;
for k = 2:(M-1);
T(j,k) = (alpha*dt)*((T(j-1,k+1)-(2*T(j-1,k))+T(j-1,k-1))/((xm/(M1))^2))+T(j-1,k);
end
end
Plotting for all 3 M values on a grid:
T(i,N) vs. xi
300
M=10
M=20
M=40
290
280
270
T(i,N)
260
250
240
230
220
210
200
0.02
0.04
0.06
0.08
0.1
xi
Comparison of 3 values of M:
In terms of convergence M=10 converged the fastest to the initial temperature 300,
followed by M=20 and M=40 for a time difference equal to 1.14s (dt=0.14s).
As M is increases from the above graph it can be seen that temperature profile becomes
smoother along the length of the beam.
