Tutorial 7

1. The velocity of an object, travelling along a straight line, was measured at various times as
follows:

Time 0 1 2 3 4 5 6 7 8 9 10 11 12
(min)

Velocity 0.0 0.6 1.7 3.4 6.3 11.18 19.09 32.12 53.60 89.02 147.41 243.69 402.43
(cm/min) 0 5 2 8 9

Estimate the distance travelled in 12 minutes using (i) Simpson’s 1/3 rule, (ii) Simpson’s 3/8 rule,
and (iii) Romberg integration, applied to Trapezoidal rule estimates with h=1, 2, and 4.

Solution:

t v Simp 1/3 Simp 3/8 Trap h=1 Trap h=2 Trap h=4
0 0 0.33 1.72 12.78
1 0.65 1.44 3.97 1.19
2 1.72 2.60 8.11
3 3.48 7.34 4.94
4 6.39 28.23 8.79 25.48 119.98
5 11.18 23.40 15.14
6 19.09 25.61 72.69
7 32.12 67.06 136.98 42.86
8 53.6 71.31 201.01 912.06
9 89.02 185.70 118.22
10 147.41 624.28 195.55 549.84
11 243.69 508.20 323.06
12 402.43
793.14 793.46 809.57 858.85 1044.82

Romberg O(h4) 793.14 796.86

Romberg O(h6) 792.89

The distance travelled in 12 minutes is: 793.14 cm using Simpson’s 1/3 rule, 793.46 cm using
3/8 rule, and 792.89 cm using Romberg integration with three trapezoidal estimates (809.57
using h=1, 858.85 using h=2, 1044.82 using h=4). Note that the Romberg, O(h4), with h=1 and
2, is same as Simpson’s 1/3. Also note that 1/3 rule is more accurate than 3/8 (assuming
Romberg to be the best estimate).

4
2. The flow rate through a circular pipe is given by . = /6 5 201 2 31, where v is the velocity at a
distance of r from the centre of pipe and r0 is the radius of the pipe. If the velocity is
4 :/;
approximated by 2 = 2 71 8 4 9 (in m/s), and the pipe radius is 12 cm, compute Q using
5
 Gauss-Legendre quadrature with 2, 3, and 4 Gauss points. Perform an error analysis for the
results, using the true value of the flow rate as 0.07389025921170497 m3/s.
Solution:

Conversion to standard domain: z=r/0.0681 => I=0.06 Iz

4 :/;
The table below shows the computation of Iz and I: @ = 201 A 2 71 8 6.:B9

2-point z w r f w.f T.V. 0.07389

-0.57735 1 0.025359 0.308044 0.308044
0.57735 1 0.094641 0.952475 0.952475 I Error (%)
Iz 1.260519 0.075631 -2.35605

3-point z w r f w.f
-0.77460 0.555556 0.013524 0.167072 0.092818
0 0.888889 0.060000 0.682900 0.607022
0.77460 0.555556 0.106476 0.979540 0.544189 I Error (%)
Iz 1.244028 0.074642 -1.01697

4-point z w r f w.f
-0.86114 0.347855 0.008332 0.103630 0.036048
-0.33998 0.652145 0.039601 0.469970 0.306489
0.33998 0.652145 0.080399 0.862339 0.56237
0.86114 0.347855 0.111668 0.958623 0.333462 I Error (%)
Iz 1.238368 0.074302 -0.55737
 Tutorial 8
1. Solve the differential equation dy/dt = −100 y + 99 e−t with the initial condition y(0)=2
using, (a) Euler’s forward (explicit) method, and (b) Euler backward (implicit)
method, to obtain the value of y at t=0.1. Use time steps of 0.01, 0.02 and 0.025. Find
the analytical solution and compare the errors for these time steps.

Solution:

The analytical solution is y= e−t+ e−100t. Formulae used:

Euler Forward: yn +1 = y n + hf (t n , yn ) Euler Backward: yn +1 = y n + hf (t n +1 , yn +1 )

Euler Explicit: Implicit Euler:

Delta t 0.01 Delta t 0.01
t y f Exact y Error t y Exact y Error
0 2 -101 2 0 0 2 2 0
0.01 0.99 -0.98507 1.35793 0.36793 0.01 1.49007 1.35793 -0.13215
0.02 0.98015 -0.97526 1.11553 0.13538 0.02 1.23024 1.11553 -0.1147
0.03 0.9704 -0.96556 1.02023 0.04984 0.03 1.09549 1.02023 -0.07526
0.04 0.96074 -0.95595 0.97911 0.01836 0.04 1.02333 0.97911 -0.04423
0.05 0.95118 -0.94644 0.95797 0.00679 0.05 0.98253 0.95797 -0.02456
0.06 0.94172 -0.93702 0.94424 0.00253 0.06 0.95744 0.94424 -0.01319
0.07 0.93235 -0.9277 0.93331 0.00096 0.07 0.94025 0.93331 -0.00695
0.08 0.92307 -0.91847 0.92345 0.00038 0.08 0.92707 0.92345 -0.00362
0.09 0.91389 -0.90933 0.91405 0.00017 0.09 0.91593 0.91405 -0.00188
0.1 0.90479 -0.90028 0.90488 9.1E-05 0.1 0.90586 0.90488 -0.00098

Delta t 0.02 Delta t 0.02

t y f Exact y Error t y Exact y Error
0 2 -101 2 0 0 2 2 0
0.02 -0.02 99.0397 1.11553 1.13553 0.02 1.3136 1.11553 -0.19806
0.04 1.96079 -100.961 0.97911 -0.98169 0.04 1.07199 0.97911 -0.09288
0.06 -0.05843 99.0777 0.94424 1.00267 0.06 0.97889 0.94424 -0.03465
0.08 1.92312 -100.924 0.92345 -0.99967 0.08 0.93555 0.92345 -0.0121
0.1 -0.09535 99.1143 0.90488 1.00024 0.1 0.90904 0.90488 -0.00416

Delta t 0.025 Delta t 0.025

t y f Exact y Error t y Exact y Error
0 2 -101 2 0 0 2 2 0
0.025 -0.525 149.056 1.05739 1.58239 0.025 1.26111 1.05739 -0.20372
0.05 3.20139 -225.967 0.95797 -2.24342 0.05 1.03297 0.95797 -0.07501
0.075 -2.4478 336.626 0.9283 3.37609 0.075 0.95118 0.9283 -0.02289
0.1 5.96786 -507.207 0.90488 -5.06298 0.1 0.91162 0.90488 -0.00673

Clearly, explicit method is unstable for step size of 0.025, and is oscillating in a finite
band for 0.02.

2. The amount of lowering of water level, s, in a well at a time t, due to pumping from
groundwater is governed by an equation of the form s=C W(u), where C is a constant
(proportional to the discharge), W is called the Well Function, and u is inversely
proportional to t. The well function is given by the equation dW(u)/du = −Exp(−u)/u.
If the value of W(1) is 0.21938393, find the value of W(0.5) using (a) Romberg
integration algorithm with accuracy O(h6) (b) Modified Euler with h= −0.5 (c) Heun’s
method with h= −0.25 and (d) Fourth-order Runge-Kutta method with h= −0.25.
Perform an error analysis using the true value of W(0.5) as 0.55977359.

Solution:
: C DE
Romberg: Use trapezoidal with h=0.125, 0.25, 0.5 to estimate /6.H 8 F
3G
 h h 
Modified Euler: y n +1 = y n + hf  t n + , y n + f (t n , y n )
 2 2 

f (t n , yn ) + f (t n +1 , yn + hf (t n , yn ))
Heun’s method: yn +1 = y n + h
2
 1 1 
k1 = f (t n , yn ); k 2 = f  t n + h, yn + hk1 
 2 2 
4th order R-K method:
 1 1 
k3 = f  t n + h, yn + hk 2 ; k 4 = f (t n + h, yn + hk3 )
 2 2 
h
yn +1 = y n + (k1 + 2k 2 + 2k3 + k 4 )
6
f(u)=dW/du=-Exp(-u)/u Romberg I=From 0.5 to 1 Modified Euler

u f(u) O(h^2) W(0.75) 0.3113538

1 -0.36788 -0.39524 h=0.5 W(0.5) 0.534295
0.875 -0.47641 -0.35507 h=0.25 Error 0.02548
0.75 -0.62982 -0.34414 h=0.125
0.625 -0.85642 Heun's
0.5 -1.21306 O(h^4)
-0.34169 h=0.5, 0.25 W(0.75)0 0.3113538
W(1.0)= 0.219384 -0.34050 h=0.25, 0.125 W(0.75) 0.3440966
W(0.5)= 0.559774
W(0.75)= 0.340332 O(h^6) W(0.5)= 0.559801 W(0.5)0 0.50155
-0.34042 Error -0.00003 W(0.5) 0.574457
Error -0.01468

4th order R-K h= -0.25

k1 -0.36788 -0.62982
k2 -0.47641 -0.85642
k3 -0.47641 -0.85642
k4 -0.62982 -1.21306

W(0.75) 0.340357 0.559880

Error -0.00011
 Tutorial 9
1. Solve the differential equation dy/dx = x2y − 2y with y(0)=1 over the interval x=0 to 0.5, using
(a) Heun’s method without iteration with h=0.25 and 0.125, (b) Heun’s method with iteration (with
h=0.25 and stopping error criterion of 1%), and (c) 4th order Runge-Kutta method with h=0.125 and
0.25. Obtain the exact value of y at x=0.5 and perform an error analysis.

Solution: f(x,y) = x2y − 2y. Exact solution is y=exp(x3/3-2x) and value at 0.5 is 0.383531573.

Heun's Without iteration

h= 0.25
xi yi f(xi,yi) y0 xi+1 f(xi+1,y0) yi+1 Error
0 1 -2 0.5 0.25 -0.968750 0.62890625
0.25 0.6289063 -1.21851 0.32428 0.5 -0.567490 0.40565681 -0.02213

h= 0.125
xi yi f(xi,yi) y0 xi+1 f(xi+1,y0) yi+1
0 1 -2 0.75 0.125000 -1.488281 0.78198242
0.125000 0.7819824 -1.55175 0.58801 0.250000 -1.139277 0.61379344
0.250000 0.6137934 -1.18922 0.46514 0.375000 -0.864870 0.48541249
0.375000 0.4854125 -0.90256 0.37259 0.500000 -0.652036 0.388250 -0.004718

Heun's With iteration

h= 0.25
xi yi f(xi,yi) y0 xi+1 f(xi+1,y0) yi+1 Error (%)
0 1 -2 0.5 0.25 -0.968750 0.62890625
0 1 -2 0.62891 0.25 -1.218506 0.59768677 -5.223385
0 1 -2 0.59769 0.25 -1.158018 0.60524774 1.2492353
0 1 -2 0.60525 0.25 -1.172667 0.60341656 -0.303467

0.25 0.6034166 -1.16912 0.31114 0.5 -0.544489 0.38921547 20.06056

0.25 0.6034166 -1.16912 0.38922 0.5 -0.681127 0.37213573 -4.589653
0.25 0.6034166 -1.16912 0.37214 0.5 -0.651238 0.37587192 0.9940069 0.007660

4th order R-K

h= 0.25
x 0 0.125 0.125 0.25
y 1.000000 0.75 0.81396 0.5962
k -2.000000 -1.48828 -1.61521 -1.15513 y(0.25)= 0.609912

x 0.25 0.375 0.375 0.5

y 0.609912 0.4621988 0.50249 0.37633 Error
k -1.181704 -0.85940 -0.93431 -0.65858 y(0.5)= 0.383757 -0.0002255
4th order R-K
h= 0.125
x 0 0.0625 0.0625 0.125
y 1.000000 0.875 0.89084 0.77773
k -2.000000 -1.74658 -1.77820 -1.54330 y(0.125)= 0.779315

x 0.125 0.1875 0.1875 0.25

y 0.779315 0.6826621 0.69548 0.6085
k -1.546454 -1.34132 -1.36651 -1.17897 y(0.25)= 0.609709

x 0.25 0.3125 0.3125 0.375

y 0.609709 0.5358772 0.546 0.47988
k -1.181311 -1.01942 -1.03867 -0.89227 y(0.375)= 0.480756

x 0.375 0.4375 0.4375 0.5

y 0.480756 0.4248866 0.43273 0.38293
k -0.893905 -0.76845 -0.78263 -0.67012 y(0.5)= 0.383544 -0.0000120

2. Solve the differential equation dy/dx = 10 sin(πx) with the initial condition y(0)=0 and step length
of 0.2 using (a) the 4th order R-K method, (b) the Milne’s method and (c) 4th order Adams method to
obtain the value of y at t=0.2, 0.4, 0.6, 0.8 and 1.0. (For the multi-step methods use the values
obtained from the R-K method for start-up)

Solution: f(x,y) = 10 sin(πx). Exact solution is (10-10 cos πx)/π and the values are:

x 0.2 0.4 0.6 0.8 1

y 0.607918 2.199467 4.166731 5.758280 6.366198

4th order R-K h= 0.2

x= 0.2 0.4 0.6 0.8 1
k1 0.0000 5.8778 9.5106 9.5106 5.8779
k2 3.0902 8.0902 10.0000 8.0902 3.0902
k3 3.0902 8.0902 10.0000 8.0902 3.0902
k4 5.8778 9.5106 9.5106 5.8779 0.0000

y= 0.6080 2.1996 4.1670 5.7586 6.3666

Adams' h= 0.2
Find y(0.8), y(1)
Milne's h= 0.2 f(x,y)
Find y(0.8), y(1) y(0.8)^0 5.6623 5.8779
f(x,y) y(0.8)^1 5.7663 5.8779
y(0.8)^0 5.6710 5.8779 y(0.8)^2 5.7663 5.8779
y(0.8)^1 5.7616 5.8779
y(0.8)^2 5.7616 5.8779

y(1)^0 6.2790 0.0000 y(1)^0 6.2759 0.0000

y(1)^1 6.3684 0.0000 y(1)^1 6.3800 0.0000
y(1)^2 6.3684 0.0000 y(1)^2 6.3800
 Tutorial 10
1. Solve the differential equation d2y/dx2 − dy/dx −2y + 2x = 3 with the boundary conditions y(0)=0
and y(0.5)=0.6967 using (a) the shooting method and (b) the direct method (use ∆x=0.25 for both).

Solution:

(a) Note that Ralston’s method has been used to solve the IVP

f1=dy1/dx=y2
f2=dy2/dx=3+2y1-2x+y2
h= 0.25
Assume y2(0)=0 Assume y2(0)=1
x y1 y2 f1 f2 x y1 y2 f1 f2
0 0 0 0 3 0 0 1 1 4
0.1875 0 0.5625 0.5625 3.1875 0.1875 0.1875 1.75 1.75 4.75
0.25 0.09375 0.78125 0.25 0.375 2.125

x y1 y2 f1 f2 x y1 y2 f1 f2
0.25 0.09375 0.78125 0.78125 3.46875 0.25 0.375 2.125 2.125 5.375
0.4375 0.24023 1.43164 1.43164 4.03711 0.4375 0.77344 3.13281 3.13281 6.80469
0.5 0.39746 1.74316 0.5 1.07422 3.70703

Therefore, y2(0)= 0.44217

x y1 y2 f1 f2
0 0 0.44217 0.44217 3.44217
0.1875 0.08291 1.08757 1.08757 3.87838
0.25 0.21811 1.37541

x y1 y2 f1 f2
0.25 0.21811 1.37541 1.37541 4.31163
0.4375 0.476 2.18384 2.18384 5.26084
0.5 0.6967 2.61152

Note: 0.21811 may be directly obtained by linear interpolation between 0.09375 and 0.375

(b) There are three nodes: at 0, 0.25, and 0.5. Only one unknown, y(0.25).

The finite difference equation at x=0.25 is

IB 8 2I: J I6 IB 8 I6
8 8 2I: J 2 A K: = 3 => 18I6 8 34I: J 14IB = 3 8 2 A K:
0.25B 0.5

With the known values of x1=0.25, y0=0 and y2=0.6967, we get y1=0.21335 (compared to 0.21811)

(If we use h=0.125: 5 nodes, 3 unknowns, the equations are

8130I: J 60IB = 3 8 2 A 0.125 = 2.75

68I: 8 130IB J 60IL = 3 8 2 A 0.25 = 2.5

68IB 8 130IL = 3 8 2 A 0.375 8 60 A 0.6967 = 839.552

and the solution is 0.07713, 0.21294, 0.41563)

2. For a plate of size Lx=2 cm and Ly=4 cm, with the boundary conditions as T=0 for y=0 and y=4;
T=100 for x=2; and insulated boundary at x=0, find the steady state temperature at the centre using
the Finite Difference Method (use ∆x=1 cm and ∆y=1 cm).

Solution: T=0
j=4

j=3

∂T/∂x=0 T=100
j=2

j=1

j=0
i=0 i=1 i=2
T=0

The discretization is shown above. We need to find T1,2. (Six unknowns are T0,1, T0,2, T0,3, T1,1, T1,2, and
T1,3)

MN O MN O
J = 0. The nodal equation, using central difference
MP N MQ N
The governing equation is
ORST,U VBOR,U WORDT,U OR,USTVBOR,UWOR,UDT
J = 0. Since
∆P N ∆Q N
approximation for the second derivative is

∆x and ∆y are equal, 8YZ,[V: 8 YZV:,[ J 4YZ,[ 8 YZ,[W: 8 YZW:,[ = 0.

Using ghost nodes (i=−1), YV:,[ = Y:,[ . The system of equation is:

 4 − 2 −1  T0,1   0 
− 1 4 0 −1   T  100
   1,1   
− 1 0 4 − 2 −1  T0, 2   0 
   =  
 −1 −1 4 0 − 1  T1, 2  100
 −1 0 4 − 2 T0,3   0 
    
 − 1 − 1 4  T1,3  100

The solution is T0,1=37.5, T1,1=50, T0,2=50, T1,2=62.5, T0,3=37.5, and T1,3=50.

3. Toxic pollutant transport in a river is governed by the following equation:

\] \] \ B] \]
J2 8 _ B J `] = 0; 0 b K b 1; ](0, ^) = ]6 ; c d = 0; ](K, 0) = K B f VP
\^ \K \K \K (:,e)
Set-up the matrix equations for the solution of the above equation in terms of Courant Number (or
CFL number) and Grid Peclet Number for general ∆x and ∆t :

2∆^
Courant or CFL No. i =
∆K

2∆K
Grid Peclet No. lC =
_

Solution:

Using forward difference for time derivative and central for the space (time is denoted by
superscript, n, and space by subscript, i), we get (for spatial variables, using a weight of μ for the
known time step, n, and a weight of 1−μ for the unknown time step, n+1)

]ZmW: 8 ]Zm ]ZW:

mW:
8 ]ZV:
mW:
]ZW:
m
8 ]ZV:
m
J 2 n(1 8 o) Jo p
∆^ 2∆K 2∆K
]ZW:
mW:
8 2]ZmW: J ]ZV:
mW:
]ZW:
m
8 2]Zm J ]ZV:
m
8 _ n(1 8 o) J o p
∆K B ∆K B
J `q(1 8 o)]ZmW: J o]Zm r = 0

Which may be written as,

i i i i i
s(1 8 o) t8 8 uv ]ZV:
mW:
J s1 J (1 8 o) t`∆^ J 2 uv ]ZmW: J s(1 8 o) t 8 uv ]ZW:
mW:
2 lC lC 2 lC
i i i i i
= so t J uv ]ZV:
m
J s1 8 o t`∆^ J 2 uv ]Zm J so t8 J uv ]ZW:
m
2 lC lC 2 lC

Or,

wZ,ZV: ]ZV:
mW:
J wZ,Z ]ZmW: J wZ,ZW:]ZW:
mW:
= xZ

The matrix form is (given c0=0 and derivative=zero at the nth node)
 a1,1 a1, 2  c1n +1   b1 
a a2, 2 a2 , 3   n +1   
 2,1  c2   b2 
 a3 , 2 a3 , 3 a3 , 4  c3n +1   b3 
  = 
 . . .  .   . 
 an −1,n − 2 an −1,n −1 an −1, n  cnn−+11  bn −1 
    
 an, n −1 + an, n +1 an, n  cnn +1   bn 