#include <stdio.h>
 
int main(void)
{
int # กำหนดค่าคงที่สำหรับการคำนวณ
k = 32  # W/m·°C
q_dot = 1.97e6  # W/m³
h = 45  # W/m²·°C
T_inf = 24  # °C
T_fixed = 50  # °C
dx = 0.08  # m
A = k / dx**2  # ค่าคงที่ของอนุพันธ์อันดับสอง
 
# สร้างเมทริกซ์ A และเวกเตอร์ b
A_matrix = np.zeros((9, 9))
b_vector = np.zeros(9)
 
# กำหนดค่าขอบ Dirichlet (T = 50°C) สำหรับโหนด 1, 4, 7
for i in [0, 3, 6]:  # โหนด 1, 4, 7
    A_matrix[i, i] = 1
    b_vector[i] = T_fixed
 
# โหนดภายใน (โหนด 5)
A_matrix[4, 4] = -4 * A
A_matrix[4, 1] = A  # บน (โหนด 2)
A_matrix[4, 3] = A  # ซ้าย (โหนด 4)
A_matrix[4, 5] = A  # ขวา (โหนด 6)
A_matrix[4, 7] = A  # ล่าง (โหนด 8)
b_vector[4] = -q_dot
 
# โหนดที่มีการพาความร้อน (โหนด 2, 3, 6, 9)
conv_nodes = [1, 2, 5, 8]
for i in conv_nodes:
    A_matrix[i, i] = -4 * A - h  # ใช้เงื่อนไขพาความร้อน
    if i - 3 >= 0:
        A_matrix[i, i - 3] = A  # บน
    if i + 3 < 9:
        A_matrix[i, i + 3] = A  # ล่าง
    if i % 3 != 0:
        A_matrix[i, i - 1] = A  # ซ้าย
    if (i + 1) % 3 != 0:
        A_matrix[i, i + 1] = A  # ขวา
    b_vector[i] = -h * T_inf - q_dot
 
# โหนดที่มีฉนวน (โหนด 8, 9) -> ไม่มีการสูญเสียพลังงานทางด้านล่าง
A_matrix[7, 7] = -3 * A  # โหนด 8 มีขอบฉนวนด้านล่าง
A_matrix[7, 4] = A  # บน (โหนด 5)
A_matrix[7, 6] = A  # ซ้าย (โหนด 7)
A_matrix[7, 8] = A  # ขวา (โหนด 9)
b_vector[7] = -q_dot
 
A_matrix[8, 8] = -3 * A - h  # โหนด 9 มีขอบฉนวนด้านล่างและพาความร้อนด้านขวา
A_matrix[8, 5] = A  # บน (โหนด 6)
A_matrix[8, 7] = A  # ซ้าย (โหนด 8)
b_vector[8] = -h * T_inf - q_dot
 
# แก้ระบบสมการ Ax = b
T_solution = np.linalg.solve(A_matrix, b_vector)
T_solution.reshape(3, 3)
 

I2luY2x1ZGUgPHN0ZGlvLmg+CgppbnQgbWFpbih2b2lkKQp7CmludCAjIOC4geC4s+C4q+C4meC4lOC4hOC5iOC4suC4hOC4h+C4l+C4teC5iOC4quC4s+C4q+C4o+C4seC4muC4geC4suC4o+C4hOC4s+C4meC4p+C4kwprID0gMzIgICMgVy9twrfCsEMKcV9kb3QgPSAxLjk3ZTYgICMgVy9twrMKaCA9IDQ1ICAjIFcvbcKywrfCsEMKVF9pbmYgPSAyNCAgIyDCsEMKVF9maXhlZCA9IDUwICAjIMKwQwpkeCA9IDAuMDggICMgbQpBID0gayAvIGR4KioyICAjIOC4hOC5iOC4suC4hOC4h+C4l+C4teC5iOC4guC4reC4h+C4reC4meC4uOC4nuC4seC4meC4mOC5jOC4reC4seC4meC4lOC4seC4muC4quC4reC4hwoKIyDguKrguKPguYnguLLguIfguYDguKHguJfguKPguLTguIHguIvguYwgQSDguYHguKXguLDguYDguKfguIHguYDguJXguK3guKPguYwgYgpBX21hdHJpeCA9IG5wLnplcm9zKCg5LCA5KSkKYl92ZWN0b3IgPSBucC56ZXJvcyg5KQoKIyDguIHguLPguKvguJnguJTguITguYjguLLguILguK3guJogRGlyaWNobGV0IChUID0gNTDCsEMpIOC4quC4s+C4q+C4o+C4seC4muC5guC4q+C4meC4lCAxLCA0LCA3CmZvciBpIGluIFswLCAzLCA2XTogICMg4LmC4Lir4LiZ4LiUIDEsIDQsIDcKICAgIEFfbWF0cml4W2ksIGldID0gMQogICAgYl92ZWN0b3JbaV0gPSBUX2ZpeGVkCgojIOC5guC4q+C4meC4lOC4oOC4suC4ouC5g+C4mSAo4LmC4Lir4LiZ4LiUIDUpCkFfbWF0cml4WzQsIDRdID0gLTQgKiBBCkFfbWF0cml4WzQsIDFdID0gQSAgIyDguJrguJkgKOC5guC4q+C4meC4lCAyKQpBX21hdHJpeFs0LCAzXSA9IEEgICMg4LiL4LmJ4Liy4LiiICjguYLguKvguJnguJQgNCkKQV9tYXRyaXhbNCwgNV0gPSBBICAjIOC4guC4p+C4siAo4LmC4Lir4LiZ4LiUIDYpCkFfbWF0cml4WzQsIDddID0gQSAgIyDguKXguYjguLLguIcgKOC5guC4q+C4meC4lCA4KQpiX3ZlY3Rvcls0XSA9IC1xX2RvdAoKIyDguYLguKvguJnguJTguJfguLXguYjguKHguLXguIHguLLguKPguJ7guLLguITguKfguLLguKHguKPguYnguK3guJkgKOC5guC4q+C4meC4lCAyLCAzLCA2LCA5KQpjb252X25vZGVzID0gWzEsIDIsIDUsIDhdCmZvciBpIGluIGNvbnZfbm9kZXM6CiAgICBBX21hdHJpeFtpLCBpXSA9IC00ICogQSAtIGggICMg4LmD4LiK4LmJ4LmA4LiH4Li34LmI4Lit4LiZ4LmE4LiC4Lie4Liy4LiE4Lin4Liy4Lih4Lij4LmJ4Lit4LiZCiAgICBpZiBpIC0gMyA+PSAwOgogICAgICAgIEFfbWF0cml4W2ksIGkgLSAzXSA9IEEgICMg4Lia4LiZCiAgICBpZiBpICsgMyA8IDk6CiAgICAgICAgQV9tYXRyaXhbaSwgaSArIDNdID0gQSAgIyDguKXguYjguLLguIcKICAgIGlmIGkgJSAzICE9IDA6CiAgICAgICAgQV9tYXRyaXhbaSwgaSAtIDFdID0gQSAgIyDguIvguYnguLLguKIKICAgIGlmIChpICsgMSkgJSAzICE9IDA6CiAgICAgICAgQV9tYXRyaXhbaSwgaSArIDFdID0gQSAgIyDguILguKfguLIKICAgIGJfdmVjdG9yW2ldID0gLWggKiBUX2luZiAtIHFfZG90CgojIOC5guC4q+C4meC4lOC4l+C4teC5iOC4oeC4teC4ieC4meC4p+C4mSAo4LmC4Lir4LiZ4LiUIDgsIDkpIC0+IOC5hOC4oeC5iOC4oeC4teC4geC4suC4o+C4quC4ueC4jeC5gOC4quC4teC4ouC4nuC4peC4seC4h+C4h+C4suC4meC4l+C4suC4h+C4lOC5ieC4suC4meC4peC5iOC4suC4hwpBX21hdHJpeFs3LCA3XSA9IC0zICogQSAgIyDguYLguKvguJnguJQgOCDguKHguLXguILguK3guJrguInguJnguKfguJnguJTguYnguLLguJnguKXguYjguLLguIcKQV9tYXRyaXhbNywgNF0gPSBBICAjIOC4muC4mSAo4LmC4Lir4LiZ4LiUIDUpCkFfbWF0cml4WzcsIDZdID0gQSAgIyDguIvguYnguLLguKIgKOC5guC4q+C4meC4lCA3KQpBX21hdHJpeFs3LCA4XSA9IEEgICMg4LiC4Lin4LiyICjguYLguKvguJnguJQgOSkKYl92ZWN0b3JbN10gPSAtcV9kb3QKCkFfbWF0cml4WzgsIDhdID0gLTMgKiBBIC0gaCAgIyDguYLguKvguJnguJQgOSDguKHguLXguILguK3guJrguInguJnguKfguJnguJTguYnguLLguJnguKXguYjguLLguIfguYHguKXguLDguJ7guLLguITguKfguLLguKHguKPguYnguK3guJnguJTguYnguLLguJnguILguKfguLIKQV9tYXRyaXhbOCwgNV0gPSBBICAjIOC4muC4mSAo4LmC4Lir4LiZ4LiUIDYpCkFfbWF0cml4WzgsIDddID0gQSAgIyDguIvguYnguLLguKIgKOC5guC4q+C4meC4lCA4KQpiX3ZlY3Rvcls4XSA9IC1oICogVF9pbmYgLSBxX2RvdAoKIyDguYHguIHguYnguKPguLDguJrguJrguKrguKHguIHguLLguKMgQXggPSBiClRfc29sdXRpb24gPSBucC5saW5hbGcuc29sdmUoQV9tYXRyaXgsIGJfdmVjdG9yKQpUX3NvbHV0aW9uLnJlc2hhcGUoMywgMykK