i s t i a r t o

import numpy as np


# difusi 1D, metode beda hingga, skema implisit

# noinspection PyShadowingNames
def algoritma_thomas(a, b, c, d):
    n = len(d)
    # matriks asli disimpan ke matriks baru, lokal di sini
    cp = c.copy()
    bp = b.copy()
    dp = d.copy()
    # forward sweep
    for i in range(1, n):
        m = a[i - 1] / bp[i - 1]
        bp[i] = bp[i] - m * cp[i - 1]
        dp[i] = dp[i] - m * dp[i - 1]
    # back substitution
    x = np.zeros(n)
    x[-1] = dp[-1] / bp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (dp[i] - cp[i] * x[i + 1]) / bp[i]

    return x


# data
panjang = 10.
k = 0.835
Delta_x = 1.
Delta_t = 0.1
jnode = int(panjang / Delta_x) - 1 # jumlah titik hitung
x = np.linspace(0, panjang, jnode + 2)  # posisi titik dari pangkal x=0 s.d. ujung x=10 cm

# variabel skalar dalam contoh kuliah adalah T (temperatur)
# di program ini, variabel temperatur T diganti phi agar tidak rancu dengan variabel waktu t
# syarat awal
phi = np.zeros(jnode)   # syarat awal, semua T = 0
# syarat batas
phi_batas_kiri = 100.   # syarat batas di pangkal (kiri) T = 100
phi_batas_kanan = 50.   # syarat batas di ujung (kanan) T = 50
# koefisien matriks sistem persamaan linear
aki = np.zeros(jnode - 1) # koefisien T_i-1
aka = np.zeros(jnode - 1) # koefisien T_i+1
ai = np.zeros(jnode)  # koefisien T_i
rhs = np.zeros(jnode) # konstanta di sisi kanan =
# titik hitung ke-2 s.d. 1 dari batas kanan: i = 1, ..., jnode-2 (karena indeks di python mulai 0)
koef_ki = k * Delta_t / Delta_x / Delta_x
koef_ka = koef_ki
for i in range(1, len(ai) - 1):
    aki[i - 1] = -1 * koef_ki
    aka[i] = -1 * koef_ka
    ai[i] = 1. + koef_ki + koef_ka
    rhs[i] = phi[i]
# titik hitung pertama: i = 0 (node bertetangga dengan batas kiri)
i = 0
aka[i] = -1 * koef_ka
ai[i] = 1. + koef_ki + koef_ka
rhs[i] = phi[i] + koef_ki * phi_batas_kiri
# titik hitung paling kanan: i = jnode-1 (node bertetangga dengan batas kanan)
i = len(ai) - 1
aki[i - 1] = -1 * koef_ki
ai[i] = 1. + koef_ki + koef_ka
rhs[i] = phi[i] + koef_ka * phi_batas_kanan
# tulis judul kolom di baris pertama
tulis_teks = " ".join(f"{jarak:9.1f}" for jarak in x)
print(f"    {tulis_teks}")
# iterasi phi[i] terhadap waktu t --> t+dt
t = 0.
tulis_teks = " ".join(f"{angka:9.4f}" for angka in phi)
print(f"{t:6.1f} {phi_batas_kiri:9.4f} {tulis_teks} {phi_batas_kanan:9.4f}")
for n in range(1, 201):
    t = n * Delta_t
    phi = algoritma_thomas(aki, ai, aka, rhs)
    if (t % 1) == 0:
        tulis_teks = " ".join(f"{angka:9.4f}" for angka in phi)
        print(f"{t:6.1f} {phi_batas_kiri:9.4f} {tulis_teks} {phi_batas_kanan:9.4f}")
    rhs[0] = phi[0] + koef_ki * phi_batas_kiri
    for i in range(1, len(rhs) - 1):
        rhs[i] = phi[i]
    rhs[len(rhs) - 1] = phi[len(rhs) - 1] + koef_ka * phi_batas_kanan
print("selesai")

import numpy as np


# TUGAS Makul SNPD
# Difusi 1D, skema implisit, metode volume hingga

# noinspection PyShadowingNames
def algoritma_thomas(a, b, c, d):
    n = len(d)
    # matriks asli disimpan ke matriks baru, lokal di sini
    cp = c.copy()
    bp = b.copy()
    dp = d.copy()
    # forward sweep
    for i in range(1, n):
        m = a[i - 1] / bp[i - 1]
        bp[i] = bp[i] - m * cp[i - 1]
        dp[i] = dp[i] - m * dp[i - 1]
    # back substitution
    x = np.zeros(n)
    x[-1] = dp[-1] / bp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (dp[i] - cp[i] * x[i + 1]) / bp[i]

    return x


# data
panjang = 10.   # cm
Gamma = 0.835   # cm^2/s
Delta_x = 1.    # cm
Delta_t = 0.1   # s
diameter = 0.6  # cm
luas = 0.25 * np.pi * diameter ** 2
Delta_V = luas * Delta_x
jmlvc = int(panjang / Delta_x)  # jumlah volume kontrol
Delta_x = panjang / jmlvc       # Delta_x dihitung ulang untuk penyesuaian dengan jumlah volume kontrol
x1 = np.arange(Delta_x/2, panjang, Delta_x)   # jarak/posisi volume kontrol
x = np.concatenate(([0.], x1, [10.]))         # jarak pangkal (x=0), volume kontrol, ujung (x=10)

# syarat awal
phi = np.zeros(jmlvc)
# syarat batas
phi_kiri = 100.  # der C
phi_kanan = 50.  # der C

# koefisien matriks tridiagonal
aW = np.zeros(jmlvc - 1)   # diagonal bawah: jumlah koefisien adalah jumlah vk-1
aE = np.zeros(jmlvc - 1)   # diagonal atas: jumlah koefisien adalah jumlah vk-1
aP = np.zeros(jmlvc)       # diagonal utama: jumlah koefisien adalah jumlah vk
rhs = np.zeros(jmlvc)      # konstanta di sisi kanan tanda kesamaan "="

# volume kontrol 2, 3, ..., jmlvc-1: i=1, 2, ..., jmlvc-2
koef_W = Gamma * luas * Delta_t / Delta_x / Delta_V
koef_E = koef_W
for i in range(1, jmlvc - 1):
    aW[i - 1] = -1 * koef_W
    aE[i] = -1 * koef_E
    aP[i] = koef_W + koef_E + 1.
    rhs[i] = phi[i]
# volume kontrol 1: i=0
koef_kiri = Gamma * luas * Delta_t / (Delta_x / 2) / Delta_V
aE[0] = -1 * koef_E
aP[0] = koef_kiri + koef_E + 1.
rhs[0] = phi[0] + koef_kiri * phi_kiri
# volume kontrol paling kanan: i=jmlvc-1
koef_kanan = Gamma * luas * Delta_t / (Delta_x / 2) / Delta_V
aW[jmlvc - 2] = -1 * koef_W
aP[jmlvc - 1] = koef_W + koef_kanan + 1.
rhs[jmlvc - 1] = phi[jmlvc - 1] + koef_kanan * phi_kanan

# iterasi phi[i] terhadap waktu t --> t+dt
t = 0.
tulis_teks = " ".join(f"{angka:9.1f}" for angka in x)   # Mengubah angka x menjadi teks
print(f"    {tulis_teks}")   # Menuliskan posisi volume kontrol di baris pertama
tulis_teks = " ".join(f"{angka:9.4f}" for angka in phi)     # Mengubah angka phi menjadi teks
print(f"{t:6.1f} {phi_kiri:9.4f} {tulis_teks} {phi_kanan:9.4f}")   # Menuliskan nilai awal phi
for n in range(1, 201):   # Menghitung phi pada waktu t+Delta_t
    t = n * Delta_t
    phi = algoritma_thomas(aW, aP, aE, rhs)   # Penyelesaian persamaan linear, matriks tridiagonal
    if (t % 1) == 0:   # Mencetak hasil hitungan tiap detik
        tulis_teks = " ".join(f"{angka:9.4f}" for angka in phi)
        print(f"{t:6.1f} {phi_kiri:9.4f} {tulis_teks} {phi_kanan:9.4f}")
    rhs[0] = phi[0] + koef_kiri * phi_kiri
    for i in range(1, jmlvc - 1):
        rhs[i] = phi[i]
    rhs[jmlvc - 1] = phi[jmlvc - 1] + koef_kanan * phi_kanan
print("selesai")