Python runs in your browser. Heavy or infinite loops may freeze your browser tab. Use "Stop" if needed; for heavy jobs, run locally.

Program 1

Laplace Equation

Laplace's equation on a square grid

$$ \frac{\partial^2 U}{\partial x^2} + \frac{\partial^2 U}{\partial y^2} = 0 $$

Dirichlet Boundary conditions are:

Left B.C.: $ U(0, y) = \frac{y}{1 + y^2}, \, \text{for} \, 0 < y < 1 $
Right B.C.: $ U(1, y) = \frac{y}{4 + y^2}, \, \text{for} \, 0 < y < 1 $
Bottom B.C.: $ U(x, 0) = 0, \, \text{for} \, 0 < x < 1 $
Top B.C.: $ U(x, 1) = \frac{1}{(1 + x^2) + 1}, \, \text{for} \, 0 < x < 1 $

Exact Analytical Solution:

$$ U_A = \frac{y}{(1 + x^2) + y^2} $$

import numpy as np
import matplotlib.pyplot as plt

#parameters
nx = 21  # Number of grid points in x-direction
ny = 21  # Number of grid points in y-direction
Lx = 1.0  # Length in x-direction
Ly = 1.0  # Length in y-direction
dx = Lx / (nx - 1)  # Spatial step size in x-direction
dy = Ly / (ny - 1)  # Spatial step size in y-direction
tol = 1e-6  # Tolerance for iterative solver

# Create grid points
x = np.linspace(0, Lx, nx)
y = np.linspace(0, Ly, ny)

# Initialize solution array
U = np.zeros((nx, ny))

# Apply boundary conditions
U[0, :] = y / (1 + y**2)         # Left boundary: U(0, y) = y / (1 + y^2)
U[-1, :] = y / (4 + y**2)        # Right boundary: U(1, y) = y / (4 + y^2)
U[:, 0] = 0.0                    # Bottom boundary: U(x, 0) = 0.0
U[:, -1] = 1 / ((1 + x)**2 + 1)  # Top boundary: U(x, 1) = 1 / [(1 + x^2) + 1]

# Iterative solver
while True:
    U_old = np.copy(U)
    for i in range(1, ny - 1):
        for j in range(1, nx - 1):
            U[i, j] = 0.25 * (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1])
    # Check convergence
    if np.max(np.abs(U - U_old)) < tol:
        break

# Exact analytical solution
UA = np.zeros((nx, ny))
for i in range(nx):
    for j in range(ny):
        UA[i, j] = y[j] / ((1 + x[i])**2 + y[j]**2)

plt.figure(figsize=(10,10))
plt.subplot(211)
# Plot numerical and exact solutions for fixed y = 0.5
y_fix = 0.5
y_index = np.abs(y - y_fix).argmin()  # Find the index closest to y_fix
plt.scatter(x, U[:, y_index], label="Numerical Solution", lw=2, ls='-.')
plt.plot(x, UA[:, y_index], label="Exact Solution", lw=2)
plt.xlabel("x")
plt.ylabel("U")
plt.legend()
plt.title(f"Solution for fixed y = {y_fix}")

plt.subplot(212)
# Plot numerical and exact solutions for fixed x = 0.5
x_fix = 0.5
x_index = np.abs(x - x_fix).argmin()  # Find the index closest to x_fix
plt.scatter(y, U[x_index, :], label="Numerical Solution", lw=2, ls='-.')
plt.plot(y, UA[x_index, :], label="Exact Solution", lw=2)
plt.xlabel("y")
plt.ylabel("U")
plt.legend()
plt.title(f"Solution for fixed x = {x_fix}")
plt.tight_layout()
plt.savefig('plot.png')
plt.show()

Output 1

Laplace's Equation

Using the iterative solver (such as the Gauss-Seidel method), we begin with Laplace's equation in two dimensions:

\[ \frac{\partial^2 U}{\partial x^2} + \frac{\partial^2 U}{\partial y^2} = 0. \]

Discretization of Laplace's Equation

The 2D Laplace equation is discretized using a finite difference method. In the code, this is done on a uniform grid, where the second derivatives with respect to $x$ and $y$ are approximated using central differences. Let's break this down.

Second derivative with respect to $x$:

\[ \frac{\partial^2 U}{\partial x^2} \approx \frac{U_{i+1,j} - 2U_{i,j} + U_{i-1,j}}{(\Delta x)^2} \] where $U_{i,j}$ is the value of $U$ at grid point $(i,j)$ and $\Delta x = \frac{L_x}{n_x - 1}$.

Second derivative with respect to $y$:

\[ \frac{\partial^2 U}{\partial y^2} \approx \frac{U_{i,j+1} - 2U_{i,j} + U_{i,j-1}}{(\Delta y)^2} \] where $\Delta y = \frac{L_y}{n_y - 1}$.

Discretized Laplace Equation

The discrete form of the Laplace equation (by adding the two central differences) becomes:

\[\small \frac{U_{i+1,j} - 2U_{i,j} + U_{i-1,j}}{(\Delta x)^2} + \frac{U_{i,j+1} - 2U_{i,j} + U_{i,j-1}}{(\Delta y)^2} = 0. \]

To simplify, we assume the grid is uniform, meaning $\Delta x = \Delta y = h$. This gives:

\[\small \frac{U_{i+1,j} - 2U_{i,j} + U_{i-1,j}}{h^2} + \frac{U_{i,j+1} - 2U_{i,j} + U_{i,j-1}}{h^2} = 0. \]

Multiplying through by $h^2$ and solving for $U_{i,j}$:

\[ U_{i,j} = \frac{1}{4} \left( U_{i+1,j} + U_{i-1,j} + U_{i,j+1} + U_{i,j-1} \right). \]

Recurrence Relation

The recurrence relation used in the iterative solver is derived from the above discretized equation. The updated value of $U_{i,j}$ is obtained by averaging the neighboring points:

\[ U_{i,j}^{(new)} = \frac{1}{4} \left( U_{i+1,j}^{(old)} + U_{i-1,j}^{(old)} + U_{i,j+1}^{(old)} + U_{i,j-1}^{(old)} \right). \]

Partial Differential Equations