Tutorial-03

Shreya Gupta (NA21B069)

LSFD for 2D Poisson Equation

Problem Definition
The 2D Poisson equation is solved on the unit domain Ω = [0, 1]2 :

∇2 u = g(x, y), (x, y) ∈ Ω

with source term, exact solution, and boundary condition:

g(x, y) = −2π 2 sin(πx) sin(πy), uexact (x, y) = 1 + x + sin(πx) sin(πy), u|∂Ω = 1 + x

1 import numpy as np
2 import matplotlib . pyplot as plt
3
4 # - - - - - - - - - - - - - - - - Problem definitions - - - - - - - - - - - - - - - -
5 g = lambda X , Y : -2.0*( np . pi **2) * np . sin ( np . pi * X ) * np . sin ( np . pi * Y )
6 uexc = lambda X , Y : 1.0 + X + np . sin ( np . pi * X ) * np . sin ( np . pi * Y )
7 bc = lambda X : 1 + X

Weight Function
A compact quartic spline weight is used:
(q
4
π
(1 − r̄2 )4 r̄ ≤ 1 r
w(r̄) = where r̄ =
0 r̄ > 1 rs

1 # - - - - - - - - - - - - - - - - Weighted LSFD pieces - - - - - - - - - - - - - - - - -

2 def compact_weight ( r_bar ) :
3 w = np . zeros_like ( r_bar )
4 mask = ( r_bar <= 1.0)
5 t = r_bar [ mask ]
6 w [ mask ] = np . sqrt (4/ np . pi ) * (1.0 - t * t ) **4
7 return w

LSFD Coefficients
For each interior node i with neighbourhood Si , the Laplacian operator coefficients α are
computed via a weighted least-squares fit of the Taylor expansion truncated to second
order.
 1 def lsfd_coeffs ( center_xy , neigh_xy , h ) :
2 dx = neigh_xy [: , 0] - center_xy [0]
3 dy = neigh_xy [: , 1] - center_xy [1]
4 r = np . sqrt ( dx * dx + dy * dy ) / ( h + 1e -15)
5 w = compact_weight ( r )
6 W = np . diag ( w )
7
8 S = np . vstack ([ dx , dy , 0.5* dx * dx , dx * dy , 0.5* dy * dy ]) . T
9 D = np . diag ([ h ** -1 , h ** -1 , h ** -2 , h ** -2 , h ** -2])
10 S_hat = S @ D
11
12 d = np . array ([0.0 , 0.0 , 1.0 , 0.0 , 1.0])
13 STWS = S_hat . T @ W @ S_hat
14 C = D @ np . linalg . solve ( STWS + 1e -12* np . eye (5) , S_hat . T @ W )
15

16 alpha = ( C . T @ d )
17 return alpha

System Assembly
1 def assemble_system ( Nx , Ny , min_neigh =20 , grow =1.4 , layers =3) :
2 xs = np . linspace (0.0 , 1.0 , Nx )
3 ys = np . linspace (0.0 , 1.0 , Ny )
4 X , Y = np . meshgrid ( xs , ys , indexing = ’ xy ’)
5 pts = np . column_stack ([ X . ravel () , Y . ravel () ]) # (N , 2)
6 N = pts . shape [0]
7
8 A = np . zeros (( N , N ) , dtype = float )
9 b = np . zeros (N , dtype = float )
10 h0 = min (( xs [1] - xs [0]) if Nx > 1 else 1.0 ,
11 ( ys [1] - ys [0]) if Ny > 1 else 1.0)
12 atol = 1e -12
13
14 is_bnd = (
15 np . isclose ( pts [: , 0] , 0.0 , atol = atol ) |
16 np . isclose ( pts [: , 0] , 1.0 , atol = atol ) |
17 np . isclose ( pts [: , 1] , 0.0 , atol = atol ) |
18 np . isclose ( pts [: , 1] , 1.0 , atol = atol )
19 )
20
21 for i in range ( N ) :
22 xi , yi = pts [ i ]
23 if is_bnd [ i ]:
24 A [i , i ] = 1.0
25 b[i] = bc ( xi )
26 continue
27
28 radius = layers * h0
29 for _ in range (8) :
30 dx = pts [: , 0] - xi
31 dy = pts [: , 1] - yi
32 mask = ( dx * dx + dy * dy ) <= radius * radius
33 idxs = np . where ( mask ) [0]
34 if idxs . size >= min_neigh :
35 break
36 radius *= grow

2
37
38 if i not in idxs :
39 idxs = np . unique ( np . append ( idxs , i ) )
40
41 neigh_xy = pts [ idxs ]
42 alpha = lsfd_coeffs ( np . array ([ xi , yi ]) , neigh_xy , h = radius
)
43 A [i , idxs ] += alpha
44 A [i , i ] += - np . sum ( alpha )
45 b[i] = g ( xi , yi )
46
47 return xs , ys , A , b

Solve and Plot

1 def solve_grid ( Nx , Ny ) :
2 xs , ys , A , b = assemble_system ( Nx , Ny )
3 u_vec = np . linalg . solve (A , b )
4 U = u_vec . reshape (( Ny , Nx ) )
5 X , Y = np . meshgrid ( xs , ys , indexing = ’ xy ’)
6 Ue = uexc (X , Y )
7 Err = np . abs ( U - Ue )
8 return xs , ys , X , Y , U , Ue , Err
9

10 def s ur f a ce _ p l ot s _ 10 x 1 0 () :
11 Nx = Ny = 10
12 xs , ys , X , Y , U , Ue , _ = solve_grid ( Nx , Ny )
13 fig = plt . figure ( figsize =(12 , 4) )
14 ax1 = fig . add_subplot (1 , 2 , 1 , projection = ’3 d ’)
15 s1 = ax1 . plot_surface (X , Y , U , cmap = ’ viridis ’)
16 ax1 . set_title ( ’ LSFD_METHOD (10 x10 ) ’)
17 ax1 . set_xlabel ( ’x ’) ; ax1 . set_ylabel ( ’y ’) ; ax1 . set_zlabel ( ’u ’)
18 fig . colorbar ( s1 , ax = ax1 , fraction =0.046 , pad =0.04)
19 ax2 = fig . add_subplot (1 , 2 , 2 , projection = ’3 d ’)
20 s2 = ax2 . plot_surface (X , Y , Ue , cmap = ’ viridis ’)
21 ax2 . set_title ( ’ Analytical (10 x10 ) ’)
22 ax2 . set_xlabel ( ’x ’) ; ax2 . set_ylabel ( ’y ’) ; ax2 . set_zlabel ( ’u ’)
23 fig . colorbar ( s2 , ax = ax2 , fraction =0.046 , pad =0.04)
24 plt . tight_layout ()
25 plt . show ()
26
27 def run_and_plot ( Nx , Ny , title_note = " " ) :
28 xs , ys , X , Y , U , Ue , Err = solve_grid ( Nx , Ny )
29 max_err = float ( np . max ( Err ) )
30 print ( f " [{ Nx } x { Ny }] max error = { max_err *100} " )
31 fig = plt . figure ( figsize =(15 , 4) )
32 ax = fig . add_subplot (1 , 3 , 1)
33 pc = ax . pcolormesh (X , Y , U , shading = ’ auto ’ , cmap = ’ viridis ’)
34 fig . colorbar ( pc , ax = ax )
35 ax . set_title ( f " Numerical u { title_note } " )
36 ax . set_xlabel ( " x " ) ; ax . set_ylabel ( " y " ) ; ax . set_aspect ( ’ equal ’ , ’
box ’)
37 ax = fig . add_subplot (1 , 3 , 2)
38 pc = ax . pcolormesh (X , Y , Ue , shading = ’ auto ’ , cmap = ’ viridis ’)
39 fig . colorbar ( pc , ax = ax )
40 ax . set_title ( " Analytical u " )

3
41 ax . set_xlabel ( " x " ) ; ax . set_ylabel ( " y " ) ; ax . set_aspect ( ’ equal ’ , ’
box ’)
42 ax = fig . add_subplot (1 , 3 , 3)
43 pc = ax . pcolormesh (X , Y , Err , shading = ’ auto ’ , cmap = ’ inferno ’)
44 fig . colorbar ( pc , ax = ax )
45 ax . set_title ( " | Error | " )
46 ax . set_xlabel ( " x " ) ; ax . set_ylabel ( " y " ) ; ax . set_aspect ( ’ equal ’ , ’
box ’)
47 plt . tight_layout ()
48 plt . show ()
49 return max_err
50

51 if __name__ == " __main__ " :

52 s ur f a ce _ p lo t s _ 10 x 1 0 ()
53 e1 = run_and_plot (10 , 10 , " (10 x10 ) " )
54 e2 = run_and_plot (20 , 20 , " (20 x20 ) " )
55 print ( f " Max error coarse : { e1 :.6 e } | fine : { e2 :.6 e } " )

4
Results
3D Surfaces — LSFD vs. Analytical (10×10 grid)
LSFD Method (10x10) Analytical (10x10)
2.4 2.4

2.6 2.2
2.4 2.2 2.4
2.2 2.2
2.0 u 2.0 2.0 u 2.0
1.8 1.8
1.6 1.6 1.8
1.4 1.8 1.4
1.2 1.2
1.0 1.6 1.0 1.6
1.0 1.0
0.8 1.4 0.8 1.4
0.0 0.2 0.6 0.0 0.2 0.6
0.4 y 0.4 y
0.4 0.2 1.2 0.4 0.2 1.2
0.6 0.6
x 0.8 0.0 x 0.8 0.0
1.0 1.0

Figure 1: 3D surface plots of the LSFD numerical solution (left) and the analytical solution
u = 1 + x + sin(πx) sin(πy) (right) on a 10 × 10 grid. The two surfaces are visually
indistinguishable at this resolution.

Output for 100 Nodes (10×10 grid)

[10x10] max error = 4.692804097756209
Numerical u (10x10) Analytical u |Error|
1.0 1.0 2.4 1.0
2.4
0.04
0.8 2.2 0.8 2.2 0.8

2.0 2.0 0.03

0.6 0.6 0.6
1.8 1.8
y

0.4 0.4 1.6 0.4 0.02

1.6

0.2 1.4 0.2 1.4 0.2 0.01

1.2 1.2
0.0 0.0 0.0
1.0 1.0 0.00
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x x

Figure 2: Pcolormesh plots for the 10 × 10 grid: numerical u (left), analytical u (centre),
and pointwise absolute error |u − uexact | (right). Peak error occurs near the domain centre
where sin(πx) sin(πy) is largest.

Increasing Nodes to 400 (20×20 grid)

[20x20] max error = 1.1248131414629103

Convergence: Max Absolute Error vs. Number of Nodes

5
 Numerical u (20x20) Analytical u |Error|
1.0 1.0 1.0
2.4 2.4 0.010
0.8 2.2 0.8 2.2 0.8
0.008
2.0 2.0
0.6 0.6 0.6
1.8 1.8 0.006
y

y
0.4 1.6 0.4 1.6 0.4
0.004
1.4 1.4
0.2 0.2 0.2
0.002
1.2 1.2
0.0 1.0 0.0 1.0 0.0 0.000
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x x

Figure 3: Same plots for the refined 20 × 20 grid. The max error drops from 4.69% to
1.12% — approximately a 4× reduction when the node count quadruples, consistent with
second-order spatial convergence.

LSFD Convergence: Max Error vs. Number of Nodes

10 1
Max absolute error

10 2

102 103
Number of nodes

Figure 4: Log-log plot of max absolute error against total node count. The nearly straight
line confirms algebraic (power-law) convergence of the LSFD scheme for this problem.

Parameter Summary

Parameter Value Description

Domain Ω [0, 1]2 Unit square
uexact 1 + x + sin(πx) sin(πy) Manufactured exact solution
g(x, y) −2π 2 sin(πx) sin(πy) Source term (RHS of Poisson eq.)
BC 1+x Applied on all four boundaries
min neigh 20 Minimum neighbours in support cloud
grow 1.4 Radius growth factor if cloud is too small
layers 3 Initial radius = layers ×h0
Grids tested 102 , 202 Coarse and fine uniform grids

Key observations:
ˆ LSFD approximates ∇2 u at each interior node via a weighted least-squares fit over
a local support cloud, avoiding a structured stencil.

6
ˆ Quadrupling the nodes (100 → 400) reduces the max error by ≈ 4× (4.69% →
1.12%), indicating second-order convergence.
ˆ The compact quartic weight ensures each node only interacts with its local neigh-
bourhood, keeping the global matrix A sparse.