<jupyter_start><jupyter_text>Half-Half Method<jupyter_code>import numpy
from scipy import linalg
from matplotlib import pyplot
%matplotlib inline
L = 1.0
nx = 51
dx = L / (nx - 1)
x = numpy.linspace(0.0, L, num=nx)

alpha = 1.22e-3 # thermal diffusivity

# initial temperature
T0 = numpy.zeros(nx)
T0[0] = 100.0
def solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc):
    T = T0.copy()
    N = len(T0) - 2
    A, b_bc = lhs_operator(N, dt, dx, alpha, cI, bc)
    for n in range(nt):
        b = rhs_vector(T, dt, dx, alpha, cE)
        T[1:-1] = linalg.solve(A, b - b_bc)
        update_boundaries(T, dx, bc)
    return T
def lhs_operator(N, dt, dx, alpha, cI, bc):
    I = numpy.identity(N)
    L, b_bc = laplacian_1d(N, dx, bc)
    A = I / dt - cI * alpha * L
    return A, b_bc
def rhs_vector(T, dt, dx, alpha, cE):
    b = T[1:-1] / dt + cE * alpha * (T[:-2] - 2 * T[1:-1] + T[2:]) / dx**2 
    return b
def update_boundaries(T, dx, bc):
    T[0] = get_bounday_value(T[1], dx, bc['left'])
    T[-1] = get_bounday_value(T[-2], dx, bc['right'])
    return 
def laplacian_1d(N, dx, bc):
    D = numpy.diag(-2.0 / dx**2 * numpy.ones(N))
    L = numpy.diag(1.0 / dx**2 * numpy.ones(N-1), k=-1)
    U = numpy.diag(1.0 / dx**2 *numpy.ones(N-1), k=+1)
    A = D + U + L
    b_bc = numpy.zeros(N)
    # Left boundary
    if bc['left']['type'] == 'Dirichlet':
        b_bc[0] = bc['left']['value'] / dx**2
    elif bc['left']['type'] == 'Neumann':
        A[0, 0] = -1.0 / dx**2
        b_bc[0] = bc['left']['value'] / dx
    # Right bounday
    if bc['right']['type'] == 'Dirichlet':
        b_bc[-1] = bc['right']['value'] / dx**2
    elif bc['right']['type'] == 'Neumann':
        A[-1, -1] = -1.0 / dx**2
        b_bc[0] = bc['right']['value'] / dx
    return A, b_bc
bc_left = {'type': 'Dirichlet', 'value': 100.0}
bc_right = {'type': 'Neumann', 'value': 0.0}
bc = {'left': bc_left, 'right': bc_right}
print(bc['left']['type'])
def get_bounday_value(T_neighbor, dx, bc):
    if bc['type'] == 'Dirichlet':
        return bc['value']
    elif bc['type'] == 'Neumann':
        return T_neighbor + bc['value'] * dx
    else:
        raise ValueError('Only Dirichlet and Neumann are supported')
    return 
sigma = 0.5
dt = sigma * dx**2 / alpha
nt = 100

cI = 0.5
cE = 0.5
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.figure(figsize=(6.0, 4.0))
pyplot.xlabel('x')
pyplot.ylabel('T')
pyplot.grid()

cI = 0.5
cE = 0.5
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.plot(x, T, label='Crank-Nicolson')

cI = 0.0
cE = 1.0
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.plot(x, T, label='Explicit')

cI = 1.0
cE = 0.0
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.plot(x, T, label='Implicit')

pyplot.legend()
pyplot.xlim(0.0, L)
pyplot.ylim(0.0, 100.0);
def analytical_temperature(x, t, alpha, L, N):
    """
    Computes and returns a truncated approximation
    of the exact temperature distribution along the rod.
    
    Parameters
    ----------
    x : numpy.ndarray
        Locations at which to calculate the temperature
        as a 1D array of floats.
    t : float
        Time.
    alpha : float
        Thermal diffusivity of the rod.
    L : float
        Length of the rod.
    N : integer
        Number of terms to use in the expansion.
    
    Returns
    -------
    T : numpy.ndarray
        The truncated analytical temperature distribution
        as a 1D array of floats.
    """
    T = 100.0 * numpy.ones_like(x)
    for n in range(1, N + 1):
        k = (2 * n - 1) * numpy.pi / (2.0 * L)
        T -= (400.0 / (2.0 * L * k) *
              numpy.sin(k * x) * numpy.exp(- alpha * k**2 * t))
    return T
pyplot.figure(figsize=(6.0, 4.0))
pyplot.xlabel('x')
pyplot.ylabel('T')
pyplot.grid()

cI = 0.5
cE = 0.5
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.plot(x, T, label='Crank-Nicolson')

cI = 0.0
cE = 1.0
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.plot(x, T, label='Explicit', linewidth=6, linestyle='--')

cI = 1.0
cE = 0.0
T = solve_1d(T0, nt, dt, dx, alpha, cI, cE, bc)
pyplot.plot(x, T, label='Implicit')

T = analytical_temperature(x, nt*dt, alpha, L, N=1000)
pyplot.plot(x, T, label='analytical_temperature')

pyplot.legend()
pyplot.xlim(0.0, L)
pyplot.ylim(0.0, 100.0);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1872808043020809, 'fraction_numerical': 0.03974748655599719, 'mean_word_length': 2.6689536878216122, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '660823', 'license': 'MIT', 'n_tokens_mistral': 1903, 'n_tokens_neox': 1806, 'n_words': 487, 'path': 'Week10 Part 1 Lesson 04_05.ipynb', 'repo_name': 'Chao8219/mae6286-chao', 'size': '4277'}