CPSC 5506 EL01: Computational Science Test 1 Solutions - February 2019

This document presents solutions to a Computational Science test (CPSC 5506 EL01) covering several key concepts. The first solution involves approximating the first and second derivatives of a function using given data points and employing difference schemes. The second solution focuses on solving the 1-D heat equation using Fourier series, deriving the solution by separating variables and applying boundary conditions and initial conditions. The third solution provides a modified pseudo-code for solving a partial differential equation (PDE) and determining the stability condition of the numerical scheme. The fourth solution presents a modified Gauss-Seidel algorithm for solving a tri-diagonal square matrix. Finally, the fifth solution involves discretizing a given PDE using finite difference schemes and deriving the resulting difference equation.

Q. 1
Solution:
Given the data we determine approximate derivatives offunction f(x) at
x=2.
From the data it can be deduced that h=0.1
Difference schemes to be employed:
f 0
(x) ≈ f(x + h) − f(x − h)
2h
f 00
(x) ≈ f (x + h) + f (x − h) − 2f (x)
h2
x f (x) f 0
(x) f 00
(x)
1.9 2.981206427x − h
2.0 3.323350970 x 3.717023853−2.981206427
2∗0.1
3.717023853+2.981206427−2∗3.323350970
0.1∗0.1
= 3.67908713 = 5.151088
2.1 3.717023853x + h
Answer:Approximately, f0
(2) = 3.67908713 and f00
(2) = 5.151088

Q. 2
Solution:
Differential equation to be is solved is the 1-d heat equation:
∂u
∂t (x, t) =∂2u
∂x2 (x, t) (1)
We begin by assuming a solution of the form:
u(x, t) = X(x)T (t)
Substituting u(x,t) in (1) yields two separate ODE’s:
X 00
(x) + λX(x) = 0 (2)
T0
(t) + λT (t) = 0 (3)
where λ is a constant.
Equation (2) is a Laplace equation with general solution:
X(x) = A cos
√ λx + B sin
√ λx
The boundary conditions: ∂u
∂x (0, t) = 0
∂u
∂x (2, t) = 0
imply:A = 0 and sin
√ λ = 0 or λ = (nπ/2)2, where n=1,2,3, . . .
∴
X(x) = bn sin (n
πx
2 ) n = 1, 2, 3, . . .
Solving equation (2) gives the solution:
T (t) = cne−(nπ/2) 2t n = 1, 2, 3 . . .
Combining them and applying the principle ofsuperposition gives the full
solution:
u(x, t) = X(x)T (t) =
∞X
n=1
Bn sin (n
πx
2 )e−(nπ/2) 2t

where, Bn = bncn. This constant can be determined from the intial condition:
∂u
∂t (x, 0) = x4 − 8x2
which implies, ∞X
1
Bn sin (n
πx
2 ) = x4 − 8x2
This is the sine Fourier series for the function:f (x) = x4 − 8x2 for which the
constant Bn is given as:
Bn = 2
2
Z 2
0
(x4 − 8x2) sin (n
πx
2 ) dx (4)
So the exact solution to the DE is:
u(x, t) =
∞X
n=1
Bn sin (n
πx
2 )e−(nπ/2) 2t
with Bn as determined in (4)

Q. 3
Solution:
Given differential equation is:
∂u
∂t (x, t) =∂2u
∂x2 (x, t)
u(0, t) = t(t − 1)
u(1, t) = sin (πt)
u(x, 0) = x(x3 − 1)
(a)
We willapproximate u(x,t) with Φi,j . Index irepresents the x coordinate
and j represents t.
The modified pseudo-code which is the solution to the above equation, is:
h = 1
N
k = 1
M
Φ0,j = jk(jk − 1) for j = 0, . . . , M # boundary condition modified
ΦN,j = sin (πjk) for j = 0, . . . , M # boundary condition modified
Φi,0 = ih(i3h3 − 1) for i = 0, . . . , N # initial condition modified
for j from 0 to M-1
{
for i from 1 to N-1
{
Φi,j+1 = Φi,j + k
h2 (Φi+1,j + Φi−1,j − 2Φi,j )
}
}

(b)
From the solution loop we have:
Φi,j+1 = Φi,j + k
h2 (Φi+1,j + Φi−1,j − 2Φi,j )
= 1 − 2
k
h2 Φi,j + k
h2 (Φi+1,j − Φi−1,j )
[k = 1
M , h = 1
N ]
= 1 − 2
N 2
M Φi,j + N 2
M (Φi+1,j − Φi−1,j )
For stability of the numerical scheme, it is required that:
1 − 2
N 2
M ≥ 0
or
0 ≤ N ≤ M
√ 2
(c)
Φi,j is an approximation for u(x, t) or more precisely u(ih, jk)

Q. 4
Solution:
The given pseudo-code of Gauss-Siedel algorithm for a regular square matrix
is modified to solve a tri-diagonal square matrix.
The main loop of the algorithm is tweaked in such a way that only non-zero
elements ofthe tri-diagonalmatrix are selected in calculations ignoring all
the zero elements.Rest of the code is unchanged and willnot be produced
here.
...
for s from 1 to N do
{ old(Xs) := Xs
error = 0
X s := bs
for k from s-1 to s-1 do
{
skip if s < 1
X s := Xs − as,kX k
}
for k from s + 1 to s+1 do
{
skip if s > N
X s := Xs − as,kX k
}
X s := Xs/as,s
error = max{error|Xs − old(Xs)|}
}
...

Q. 5
Solution:
Given:
∂2u
∂x2 (x, y) +∂2u
∂y2 (x, y) = 0 (1)
The domain xy is defined as:
0 ≤ x ≤ 1;0 ≤ y ≤ 2
(a)
Index i represents x and j represents y.For M sub-intervals over each coor-
dinate length we define:
h = i
M & k = j
M
therefore,
h
k = i/M
j/M = 1
2
Let Φi,j approximate u(x,y).Discretize (1) using difference scheme:
∂2u
∂x2 (x, y) =Φi+1,j + Φi−1,j − 2Φi,j
h2
∂2u
∂y2 (x, y) =Φi,j+1 + Φi,j−1 − 2Φi,j
k2
Substituting in (1):
Φi+1,j + Φi−1,j − 2Φi,j
h2 + Φi,j+1 + Φi,j−1 − 2Φi,j
k2 = 0
Φi+1,j + Φi−1,j − 2Φi,j + h2
k2 (Φi,j+1 + Φi,j−1 − 2Φi,j ) = 0
Φi+1,j + Φi−1,j − 2Φi,j + 1
4(Φi,j+1 + Φi,j−1 − 2Φi,j ) = 0

therefore,
Φi+1,j + Φi−1,j + 0.25Φi,j+1 + 0.25Φi,j−1 − 2.5Φi,j = 0
is the required answer.
(b)
Φi,j is an approximation for u(x,y) or more precisely u(ih, jk)