แม่นก กกลูกนก (๑๒) MO Memoir : Wednesday 26 February 2563

และแล้วในที่สุด พวกเขาก็จากไป

ตอนแรกคิดว่าจะมีลูกนก ๒ ตัว แต่ท้ายสุดก็มีเพียงแค่ตัวเดียว

เช้าวันจันทร์ ก่อนออกไปทำงาน เห็นมาอยู่พร้อมหน้ากันทั้งพ่อแม่ลูก

พอตอนเย็นกลับมา ก็เห็นไปเกาะอยู่ที่สายไฟ และคงเกาะที่นั่นทั้งคืน

พอมาเช้าวันอังคาร เห็นลูกนกเกาะอยู่ตัวเดียว หยิบกล้องมาถ่ายรูปได้ ถ่ายได้เพียงรูปเดียว เขาก็บินจากไป

ค่ำวันวาน ก็เลยถือโอกาสทำความสะอาดโคมไฟ เขี่ยเอาขี้นกและเศษซากรังออกมา

เมื่อถึงเวลาที่เหมาะสม พวกเขาก็คงจะกลับมากันใหม่เอง

รูปที่ ๑ เป็นวันแรกที่เห็นพร้อมหน้าพร้อมตา

ถนนคู่ขนานทางรถไฟก็มาแล้ว ทางด่วนคร่อมทางรถไฟก็มาแล้ว รถไฟคู่ขนานทางด่วนก็รออยู่ว่าจะเริ่มเมื่อใด รถไฟลอยฟ้าที่ถนนหน้าปากซอยก็เริ่มทดลองวิ่ง พื้นที่สวนเดิมที่เคยเป็นที่ตาบอดก็ออกถนนใหญ่ได้โดยง่าย คอนโดสูงหลายหลังเริ่มปรากฏ เมื่อกลับมาทำงานใหม่ ๆ เมื่อปี ๒๕๓๗ นั้น เคยมีผู้ใหญ่ในที่ทำงานท่านหนึ่งถามว่าบ้านอยู่ไหน พอบอกว่าอยู่บางพลัด เขาก็บอกว่าอยู่ตรงนั้นได้อย่างไร หาความเจริญใด ๆ ไม่ได้เลย ซึ่งจะว่าไปมันก็เป็นจริงอย่างเขาว่า เพราะแต่ก่อนแถวนั้นไม่ค่อยมีรถผ่าน เพิ่งจะมีมากก็ตอนที่มีสะพานพระราม ๗ แล้ว บ้านเรือนริมถนนก็เป็นตึกแถวสองชั้นเล็ก ๆ อยู่ต่ำกว่าระดับถนนใหญ่อีก เรียกว่าสร้างกันก่อนที่จะมีการทำถนนใหม่ ซอยต่าง ๆ ก็เป็นซอยเล็ก ๆ ปากซอยกว้างอย่างมากก็เพียงแค่ให้รถเก๋งวิ่งสวนกันได้แค่นั้นเอง

 

รูปที่ ๒ อันที่จริงฟ้ามันมืดแล้ว แต่ถ่ายด้วยกล้องดิจิตอลเปิดแฟลช ก็เลยได้ภาพมาสว่าง

   

รูปที่ ๓ ถ่ายอีกมุมหนึ่ง

รูปที่ ๔ ภาพสุดท้ายของครอบครัวนี้ที่ถ่ายเอาไว้ได้

ช่วงหน้าหนาวที่ผ่านมา จากหน้าต่างห้องนอนลูกยังพอมองเห็นหิ่งห้อยบินเล่นในสวนข้างบ้านตอนหัวค่ำอยู่ แต่ช่วงนี้ไม่เห็นแล้ว ได้เห็นพวกมันบินเล่นอยู่ข้างบ้าน (บางทีก็แวะเข้ามาในบ้านด้วย) มากว่า ๔๐ ปี ส่วนมันจะอยู่ได้อีกนานแค่ไหนก็คงขึ้นอยู่กับว่าพื้นที่สวนตรงนั้นจะมีการเปลี่ยนแปลงอย่างไรในอนาคต

  

ส่วนกระรอกก็ยังคงวิ่งไล่จับกันตอนเช้าเหมือนเดิม

"ให้" กับ "ห้าม" ความหมายมันตรงกันข้าม MO Memoir : Monday 24 February 2563

ลองดูหัวข้อข่าวและเนื้อข่าวในรูปข้างล่างก่อนนะครับ

จากพาดหัวข่าว ถ้าว่ากันตามหลักภาษาไทย ศธ. (กระทรวงศึกษาธิการ)" คือประธาน ที่ "ออกหนังสือเวียน" โดยเนื้อหาในหนังสือเวียนคือการ "xxx" นักเรียนและครูที่เดินทางไปยัง ๖ ประเทศเสี่ยงโควิด 19 "หยุดเรียนและหยุดปฏิบัติงาน ๑๔ วัน"

  

ถ้าคำใน "xxx" คือ "ห้าม" ก็แสดงว่ากลับมาแล้วห้ามหยุด ต้องมาทำงาน ต้องมาเรียนใช่ไหมครับ

ถ้าคำใน "xxx" คือ "ให้" ก็แสดงว่ากลับมาแล้วห้ามมาทำงาน ห้ามมาเรียน ต้องหยุดใช่ไหมครับ

แล้วจะสรุปว่าอย่างไรดี :) :) :)

  

จะว่าไปประโยคในเนื้อข่าวก็อาจตีความได้อีกแบบหนึ่ง ตรงข้อความที่ว่า "...และหยุดปฏิบัติงานเฝ้าระวังอาการเป็นเวลา ๑๔ วัน ...." ทั้ง ๆ ที่สิ่งที่ควรจะเป็นมากกว่าคือ "...และหยุดปฏิบัติงาน และให้เฝ้าระวังอาการเป็นเวลา ๑๔ วัน ...." คือควรจะมีการเว้นวรรคและเพิ่มคำ เพื่อแสดงให้เห็นว่าคำว่าการหยุดนั้นมันสิ้นสุดเฉพาะ "ปฏิบัติงาน"

 

มีการเขียนแบบหนึ่งที่โดยส่วนตัวแล้วไม่ค่อยพบเจอเท่าใดนัก เว้นแต่จะเป็นในคู่มือการทำงาน คือแทนที่จะเขียนเป็นย่อหน้ายาว ๆ กลับใช้การแยกเรื่องออกเป็นข้อ ๆ แทน การเขียนแบบนี้จะช่วยลดโอกาสตีความเป็นอย่างอื่น

  

เช่นในกรณีนี้ถ้าจะเขียนใหม่ว่าเป็น

".... ขอให้นักเรียนและครูที่เดินทางไปในกลุ่มประเทศดังกล่าว

๑. หยุดเรียนและหยุดปฏิบัติงาน และ

๒. เฝ้าระวังอาการ

เป็นเวลาทั้งสิ้น ๑๔ วัน ...."

ก็น่าจะให้ความหมายที่ชัดเจนกว่า และยากที่จะตีความเป็นอย่างอื่น

วันนี้ก็ไม่มีอะไรมาก แค่เห็นประโยคที่ข่าวนี้ใช้มันน่าสนใจดี ก็เลยเอามาเขียนบันทึกไว้ เพื่อไว้เป็นตัวอย่างสอนการเขียนรายงานให้กับนิสิต :) :) :)

แม่นก กกลูกนก (๑๑) MO Memoir : Sunday 16 February 2563

บันทึกไว้หน่อยหลังจากผ่านไป ๓ สัปดาห์ว่าคราวนี้ นกเขาครอบครัวนี้มีสมาชิกเพิ่มมาอีก ๒ ตัว

จะเรียกว่ามุมนี้เป็นมุมที่เหมาะสมก็ได้ เพราะไม่มีตัวอะไรมากวนได้ เพราะวันนี้ก็เห็นคราบงูที่มาลอกเอาไว้ที่หลังคาโรงรถฝั่งตรงข้ามรังนก (รูปสุดท้าย) นี่ถ้าเป็นการทำรังบนต้นไม้ ก็คงจะเสร็จงูตัวนี้ไปแล้ว

  

รูปที่ ๑ เพิ่งจะได้ภาพลูกนกเป็นครั้งแรก เนื่องจากแม่นกออกไปหาอาหาร

นับจากที่เห็นมากกไข่ นี่ก็ผ่านไปสามสัปดาห์แล้ว แสดงว่านับจากแม่นกออกไข่ จนถึงเวลาที่ลูกนกโผบินจากรังได้ ก็กินเวลาไม่น้อยกว่าหนึ่งเดือน

  

รูปที่ ๒ อีกมุมหนึ่ง ทั้งรังเต็มไปด้วยขี้นก รอทิ้งรังเมื่อใดก็จะได้ขึ้นไปเก็บกวาด

  

รูปที่ ๓ แม่นกกลับมาป้อนอาหารให้ลูกนก 

  

รูปที่ ๔ ภาพอีกมุมหนึ่ง

   

รูปที่ ๕ หลังคาโรงรถฝั่งตรงข้ามรังนก เพิ่งจะเห็นตอนขึ้นไปถ่ายรูปลูกนก วัดความยาวคราบที่มันลอกออกมาได้ ๘๐ เซนติเมตร เดาว่าน่าจะเป็นเขียวพระอินทร์

การคำนวณเชิงตัวเลข (๒๘) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย ODE solvers ของ GNU Octave ตอนที่ ๓ MO Memoir : Saturday 15 February 2563

การเกิดปฏิกิริยาใน fixed-bed catalystic reactor นั้น ถ้า heat of reaction มีค่าต่ำ ก็จะประมาณได้ว่า reactor ตัวนั้นทำงานแบบอุณหภูมิคงที่ (isothermal) และถ้า heat of reaction มีค่าไม่สูงมาก ก็จะออกแบบให้ reactor นั้นทำงานแบบ adiabatic (ไม่มีการระบายความร้อนออกหรือให้ความร้อนแก่ reactor โดยตรง) กล่าวคือปล่อยให้ปฏิกิริยาดำเนินไปข้างหน้าเรื่อย ๆ จนได้ค่า conversion ที่ต้องการ แต่ถ้าหากว่าอุณหภูมิของระบบนั้นสูง/ต่ำเกินไปก่อนจะได้ค่า conversion ที่ต้องการ ก็ต้องมีการนำเอา reactor หลายตัวมาต่ออนุกรมกัน โดยมีเครื่องแลกเปลี่ยนความร้อน (ที่ระบาย/ป้อนความร้อน) อยู่ระหว่าง reactor แต่ละตัว

  

ในกรณีของการทำงานแบบ adiabatic นั้น ปฏิกิริยาดูดความร้อน (endothermic reaction) มีแนวโน้มที่จะหยุดตัวเอง เพราะเมื่อปฏิกิริยาดำเนินไปข้างหน้าเรื่อย ๆ ความเข้มข้นของสารตั้งต้นจะลดต่ำลง และอุณหภูมิของระบบก็จะลดต่ำลง สองปรากฏการณ์นี้ส่งผลให้อัตราการเกิดปฏิกิริยาลดต่ำลง 

    

ในกรณีของปฏิกิริยาคายความร้อน (exothermic reaction) นั้นเมื่อปฏิกิริยาดำเนินไปข้างหน้า ความเข้มข้นสารตั้งต้นที่ลดต่ำลงจะดึงให้อัตราการเกิดปฏิกิริยาลดลง แต่อุณหภูมิที่สูงขึ้นจากความร้อนที่คายออกมาจะทำให้ปฏิกิริยาเกิดเร็วขึ้น สองปรากฏการณ์นี้ให้ผลที่ขัดแย้งกัน ถ้าหากการทำปฏิกิริยานั้นเกิดขึ้นในเครื่องปฏิกรณ์แบบ adiabatic (ไม่มีการระบายความร้อนออก) ผลที่เกิดขึ้นจะเกิดในทิศทางเดียวคืออุณหภูมิของระบบจะเพิ่มสูงขึ้นเรื่อย ๆ จนอาจเหนือการควบคุม (ที่เรียกว่าปฏิกิริยาเกิดการ runaway)

   

เพื่อที่จะ "ลด" โอกาสเกิดเหตุการณ์ดังกล่าว ในกรณีของปฏิกิริยาคายความร้อนสูงจึงมักดำเนินการในรูปแบบที่เรียกว่า non-isothermal non-adiabatic คือมีการระบายความร้อนออกจาก reactor โดยตรง ดังนั้นอุณหภูมิใน reactor จึงขึ้นกับอัตราการคายความร้อนจากปฏิกิริยาและอัตราการระบายความร้อนออก เนื่องจากอัตราการระบายความร้อนออกนั้นขึ้นอยู่กับผลต่างระหว่างอุณหภูมิภายใน reactor กับอุณหภูมิแหล่งรับความร้อน ดังนั้นในช่วงที่ปฏิกิริยาเพิ่งเริ่มเกิด (ด้านขาเข้าของ reactor) อุณหภูมิใน reactor (ที่เพิ่มขึ้นจากความร้อนที่ปฏิกิริยาคายออกมา) ยังไม่สูง ผลต่างอุณหภูมิจึงยังไม่สูง อัตราการระบายความร้อนออกจึงต่ำ ส่งผลให้อุณหภูมิใน reactor เพิ่มสูงขึ้นเรื่อย ๆ เมื่อปฏิกิริยาดำเนินไปข้างหน้า แต่เมื่ออุณหภูมิใน reactor สูงถึงระดับหนึ่งจนทำให้ผลต่างอุณหภูมิสูงมากพอ อุณหภูมิใน reactor จะไม่เพิ่มขึ้นและจะเริ่มลดต่ำลง (ผลจากความเข้มข้นสารตั้งต้นที่ลดต่ำลง)

  

ตัวอย่างของปฏิกิริยาประเภทนี้ได้แก่ปฏิกิริยา partial oxidation ไฮโดรคาร์บอนไปเป็นสารประกอบ oxygenate ที่ใช้อากาศเป็นสารออกซิไดซ์ เนื่องจากออกซิเจนในอากาศมีความว่องไวต่ำในการออกซิไดซ์สารตั้งต้น การเริ่มการเกิดปฏิกิริยาจึงต้องการอุณหภูมิที่สูงมากพอ (ไม่เช่นนั้นปฏิกิริยาจะไม่เกิด หรือเกิดช้ามากจนไม่ได้ conversion ที่ต้องการ) อุณหภูมิของแหล่งรับความร้อนจึงต้องสูงตามไปด้วยเพื่อทำให้สารตั้นต้นที่ป้อนเข้ามานั้นมีอุณหภูมิสูงพอที่จะเกิดปฏิกิริยาได้ แต่ในขณะเดียวกันต้องไม่สูงจนทำให้การดึงความร้อนออกจาก reactor ต่ำเกินไป แหล่งรับความร้อนจากปฏิกิริยาที่ใช้กันจึงเป็นพวก molten salt หรือเกลือหลอมเหลวที่เกิดจากสารผสมในสัดส่วนที่เหมาะสมระหว่างเกลือบางชนิด (ทำให้มีจุดหลอมเหลวต่ำเพื่อง่ายต่อการหลอมให้เป็นของเหลว) ตัวอย่างอุณหภูมิของ molten salt นี้อาจอยู่ที่ระดับประมาณ 250-450ºC ขึ้นอยู่กับปฏิกิริยา สภาวะการทำปฏิกิริยา ความว่องไวของตัวเร่งปฏิกิริยา ความเข้มข้นสารตั้งต้น ฯลฯ

ที่ใช้คำว่า "ลด" ในย่อหน้าข้างบนนั่นก็เพราะว่าแม้ว่าจะมีการดึงความร้อนออกจาก reactor โดยตรงมันก็ยังมีโอกาสที่จะเกิดการ runaway อยู่ถ้าหากว่ามีความร้อนสะสมใน reactor มากเกินไป สิ่งที่ทำให้ความร้อนสะสมใน reactor ได้มีทั้ง ความเข้มข้นสารตั้งต้นที่สูง อัตราการไหลที่ต่ำ อุณหภูมิของแหล่งรับความร้อนที่สูงเกินไป อุณหภูมิขาเข้าที่สูงเกินไป และความว่องไวของตัวเร่งปฏิกิริยาที่สูงเกินไป ดังนั้นช่วงการทำงานของ reactor เหล่านี้จึงต้องอยู่ในช่วงที่เหมาะสมที่สามารถทำให้ได้ค่า conversion ที่สูงมากพอโดยไม่เกิดการ runaway

  

การคาดการณ์การเกิด runaway ทำได้ด้วยการแก้สมการดุลมวลสารและพลังงานของระบบ และทดลองปรับเปลี่ยนค่าพารามิเตอร์ต่าง ๆ เพื่อตรวจสอบดูว่าจะเกิดอะไรขึ้น รูปแบบง่าย ๆ รูปแบบหนึ่งของสมการดุลมวลสารและพลังงานแสดงไว้ข้างล่าง (เรียกว่าแบบจำลอง 1 มิติหรือ 1-dimensional model เพราะสมมุติให้การเปลี่ยนแปลงเกิดขึ้นตามความยาวของ reactor เท่านั้น)

  

ในที่นี้พจน์ C คือความเข้มข้น, T คืออุณหภูมิของ reactor, Tc คืออุณหภูมิของแหล่งรับความร้อน, z คือระยะทางตามความยาวของ reactor, exp(a - b/T) คือ rate constant ที่เปลี่ยนตามอุณหภูมิ, H เป็นผลรวมของความร้อนที่เกิด และ U คือค่าสัมประสิทธิ์การถ่ายเทความร้อนระหว่างภายใน reactor กับแหล่งรับความร้อน สมการดุลมวลสารนั้นประกอบด้วยพจน์ของอัตราการหายไปของสารตั้งต้น ในขณะที่สมการดุลความร้อนนั้นประกอบด้วยความร้อนที่เกิดจากปฏิกิริยาและความร้อนที่มีการดึงออก ระบบสมการทั้งสองเป็นระบบสมการอนุพันธ์สามัญที่มีพจน์ที่เป็น non-linear ร่วมอยู่

   

ในที่นี้การหาคำตอบระบบ 2 สมการข้างต้นจะใช้คำสั่ง lsode, ode45 และ ode15s ของ GNU Octave เปรียบเทียบกันโดยจะหาคำตอบในช่วง z = 0.0 - 0.2 คู่มือของโปรแกรมนั้นไม่ได้ให้รายละเอียดว่าคำสั่งดังกล่าวหาคำตอบด้วยวิธีใด คำสั่ง lsode นั้นให้กำหนดช่วงที่ต้องการหาคำตอบและระบุด้วยว่าต้องการซอยช่วงดังกล่าวออกเป็นช่วงย่อยที่ละเอียดมากน้อยเท่าใด (ยิ่งละเอียดมากก็ยิ่งเข้าใกล้คำตอบที่ถูกต้องมากขึ้น) ในที่นี้ได้ทดลองเลือกแบ่งย่อยเป็น 200 ช่วง (ในโปรแกรมต้องใส่ 201 จุด) ค่าพารามิเตอร์ต่าง ๆ ที่ใช้คือ a = 20, b = 14000, H = 10000, T ขาเข้า = 350ºC, C ขาเข้า = 1.0, และ U = 500 ตัวโปรแกรมที่เขียนขึ้นแสดงไว้ในรูปที่ ๑ ส่วนผลการคำนวณที่ได้ที่ใช้ค่า Tc = 414.0, 414.7 และ 414.8ºC แสดงไว้ในรูปที่ ๒ - ๖ โดยรูปที่ ๒ และ ๓ นั้นได้มาจากคำสั่ง lsode รูปที่ ๔ และ ๕ ได้มาจากคำสั่ง ode45 และรูปที่ ๖ ได้มาจากคำสั่ง ode15s

   

ในกรณีของคำสั่ง lsode และ ode45 นั้นพบว่าสองตัวนี้ให้คำตอบที่เหมือนกันคือการทำงานของ reactor จะเป็นถ้าอุณหภูมิของแหล่งรับความร้อนนั้นไม่เกิน 417.7ºC แต่เมื่อเพิ่มอุณหภูมิของแหล่งรับความร้อนเป็น 417.8ºC พบว่าปฏิกิริยาจะเกิดการ runaway ซึ่งเห็นได้จากผลการคำนวณนั้นให้ค่าอุณหภูมิที่เพิ่มสูงขึ้นกระทันหันโดยให้ค่าอุณหภูมิสูงสุดอยู่ที่หลายพันองศาเซลเซียส การแกว่งของคำตอบที่เกิดขึ้นกับกรณีของ lsode นั้นเกิดจากการที่จำนวนจุดคำนวณนั้นต่ำเกินไป (ใช้เพียงแค่ 200 จุด) แต่ที่ไม่พบในกรณีของ ode45 นั้นเป็นเพราะว่า ode45 ใช้การปรับระยะ step size การคำนวณโดยอัตโนมัติเพื่อให้สอดรับกับอัตราการเปลี่ยนแปลงของคำตอบ แต่สิ่งที่เห็นคือระยะห่างของจุดบริเวณนั้นแคบมาก (อยู่ที่ระดับ 10^-7-10^-6) จำนวนจุดการคำนวณตลอดช่วงอยู่ที่เกือบ 35,000 จุด ทำให้การคำนวณใช้เวลานานมากอย่างเห็นได้ชัด

  

แต่พอเปลี่ยนมาใช้ ode15s ปรากฏว่าคำตอบที่ได้นั้นเปลี่ยนไปเป็นคนละเรื่องเลย คือคำสั่งนี้ไม่สามารถพบตำแหน่งที่ปฏิกิริยาเกิดการ runaway ได้ (ตำแหน่งที่อุณหภูมิและความเข้มข้นเปลี่ยนแปลงอย่างรวดเร็ว) คำตอบที่ได้นั้นเสมือนกับว่าปฏิกิริยาแทบจะไม่เกิดแม้ว่าจะใช้อุณหภูมิ Tc สูงกว่า 415.0ºC ก็ตาม

   

รูปที่ ๑ ชุดคำสั่งที่ใช้ในการคำนวณ ข้างบนเป็นของ lsode ส่วนข้างล่างเป็นของ ode45s

ความแตกต่างของตัวอย่างนี้กับตัวอย่างในตอนที่ ๑ และ ๒ อยู่ตรงที่ในตอนที่ ๑ และ ๒ นั้นการเปลี่ยนแปลงอย่างรวดเร็วหรืออย่างกระทันหันเกิดขึ้นจากจุดเริ่มต้นการคำนวณ แต่ในกรณีนี้เกิดขึ้นระหว่างเส้นการคำนวณไปยังจุดปลายทาง ในกรณีของวิธี multi-step ที่มีการใช้จุดข้อมูลย้อนหลังในการคำนวณ ประกอบกับการเลือกระยะ step size อัตโนมัติ จึงเป็นไปได้ว่าการเปลี่ยนแปลงอย่างช้า ๆ ที่เกิดขึ้นอย่างต่อเนื่องมาก่อนหน้า ทำให้ไม่มีการคาดการณ์ว่าจะมีการเปลี่ยนแปลงอย่างกระทันหันเกิดขึ้นต่อหน้า ทำให้การคำนวณเดินข้ามตำแหน่งดังกล่าวไป 

   

ตัวอย่างนี้น่าจะเป็นอีกตัวอย่างหนึ่งที่แสดงให้เห็นความสำคัญของการทดลองวิธีการที่จะเลือกใช้แก้ปัญหาว่า ถ้าเปลี่ยนวิธีการ จะยังทำให้คำตอบที่ได้นั้นออกมาเหมือนกันอยู่หรือไม่

   

รูปที่ ๒ กราฟการเปลี่ยนแปลองอุณหภูมิภายใน reactor ที่ได้จากคำสั่ง lsode ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC โดยแบ่งช่วงการคำนวณออกเป็น 200 ช่วงย่อย

   

รูปที่ ๓ กราฟการเปลี่ยนแปลงความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง lsode ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC โดยแบ่งช่วงการคำนวณออกเป็น 200 ช่วงย่อย

   

รูปที่ ๔ กราฟการเปลี่ยนแปลองอุณหภูมิภายใน reactor ที่ได้จากคำสั่ง ode45 ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC

   

รูปที่ ๕ กราฟการเปลี่ยนแปลงความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง ode45 ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC

   

รูปที่ ๖ กราฟการเปลี่ยนแปลงอุณหภูมิและความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง ode15s ที่ค่า Tc = 415.0ºC

การคำนวณเชิงตัวเลข (๒๗) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย ODE solvers ของ GNU Octave ตอนที่ ๒ MO Memoir : Thursday 13 February 2563

กลับมาต่อกันจากเมื่อวันอังคาร ทบทวนนิดนึงว่าสมการที่นำมาใช้ทดสอบคือสมการอัตราการเกิดปฏิกิริยาของสาร A เมื่อเวลาใด ๆ ที่มีรูปแบบ โดยที่เมื่อ t = 0, A = 1 สมการนี้มีคำตอบคือ เมื่อ A₀ คือค่า A ที่เวลา t = 0 แต่คราวนี้เพิ่มค่า k เป็น 1000 หรือ ๑๐ เท่าของตอนที่แล้ว
   

การคำนวณยังคงใช้ด้วย solver ของ GNU octave ๔ ตัวด้วยกันคือ ode45, ode23, ode15s และ ode15i และเทคนิค implicit Euler ที่เขียนขึ้นมาเองที่ใช้ Δt = 0.001 คงที่ตลอดการคำนวณ และเนื่องด้วยคำตอบที่ได้นั้นมีการเปลี่ยนแปลงอย่างรวดเร็ว จึงได้เขียนกราฟถึงแค่ t = 0.01 ผลการคำนวณที่ได้แสดงไว้ในรูปที่ ๑ - ๕

   

รูปที่ ๑ คำตอบที่ได้จากการคำนวณด้วย solver ode45 เมื่อใช้ k = 1000

จะเห็นได้ว่าในกรณีนี้ ode23 ยังมีปัญหาเรื่องให้ค่าความเข้มข้นออกมาติดลบเกิดขึ้นเร็ว (แค่ประมาณ t = 0.002 ในขณะที่ ode45 และ ode15s แต่ถ้าพิจารณาความใกล้เคียงคำตอบที่ถูกต้องจะเห็นว่า ode15s ให้ค่าที่ดีกว่า ในกรณีของ ode15i นั้นก็มีปัญหาเรื่องค่าความเข้มข้นที่คำนวณได้นั้นออกมาติดลบ แต่ก็ไปเกิดขึ้นแถว ๆ t = 0.017 ส่วนในกรณีของวิธี Implicit Euler นั้นที่เห็นค่าที่คำนวณได้กับคำตอบที่ถูกต้องนั้นห่างกันมากเป็นเพราะการใช้ Δt ที่กว้าง ซึ่งถ้าลดขนาดของ Δt ลงก็จะได้คำตอบที่ใกล้เคียงกับคำตอบที่ถูกต้องมากขึ้น
    

รูปที่ ๒ คำตอบที่ได้จากการคำนวณด้วย solver ode23 เมื่อใช้ k = 1000

  

รูปที่ ๓ คำตอบที่ได้จากการคำนวณด้วย solver ode15s เมื่อใช้ k = 1000

ในทางทฤษฏีนั้นไม่ว่าจะใช้เทคนิคใด ๆ ในการแก้ปัญหา ถ้าใช้ step size ที่มีขนาดเล็กพอก็จะได้คำตอบที่ถูกต้องเหมือนกันหมด เว้นแต่ว่าจะพบกับโจทย์ที่มีปัญหาเรื่องเสถียรภาพสูงมากจนไม่สามารถใช้วิธีตระกูล explicit ได้ ทำให้ต้องใช้เทคนิคตระกูล implicit เท่านั้น
  

ตัวอย่างที่ยกมานี้ต้องการแสดงให้เห็นว่าในการเลือกใช้โปรแกรมสำเร็จรูปที่ชุดคำสั่งให้เลือกใช้นั้น ควรต้องทดลองหาคำตอบด้วยคำสั่งตอบด้วยชุดคำสั่งต่าง ๆ เพื่อดูว่าผลการคำนวณที่ได้นั้นขึ้นอยู่กับคำสั่งที่เลือกใช้หรือไม่ ถ้าพบว่าคำตอบที่ได้นั้นแตกต่างกันก็ต้องมาพิจารณาว่าอันไหนคือคำตอบที่ถูกและอันไหนคือที่ผิดเพื่อจะได้เลือกใช้ชุดคำสั่งที่เหมาะสมทำงาน

  

รูปที่ ๔ คำตอบที่ได้จากการคำนวณด้วย solver ode15i เมื่อใช้ k = 1000

   

รูปที่ ๕ คำตอบที่ได้จากการคำนวณด้วยเทคนิค Implicit Euler เมื่อใช้ k = 1000

แนวทางหัวข้อการทำวิทยานิพนธ์นิสิตรหัส ๖๑ (ตอนที่ ๘) MO Memoir : Wednesday 12 February 2563

เอกสารฉบับนี้แจกจ่ายเป็นการภายใน ไม่นำเนื้อหาลง blog

การคำนวณเชิงตัวเลข (๒๖) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย ODE solvers ของ GNU Octave ตอนที่ ๑ MO Memoir : Tuesday 11 February 2563

คู่มือการใช้งาน GNU Octave 5.1.0 นั้นกล่าวว่ามีฟังก์ชัน (ที่เรียกว่า solver) สำหรับแก้ปัญหาระบบสมการอนุพันธ์สามัญ (Ordinarty Differential Equations หรือที่ย่อกันว่า ODE) ที่ "compatible" กับของ MATLAB อยู่ด้วยกัน ๔ ตัวคือ ode45, ode23, ode15s และ ode15i (ซึ่งต่างเป็นวิธีการที่ใช้กับการแก้ไขปัญหาค่าเริ่มต้น (initial value problem) โดยวันนี้เราจะลองมาดูกันว่าผลการทำงานของแต่ละวิธีมันเป็นอย่างไร (ด้วยการทดสอบกับโจทย์ง่าย ๆ) แต่ก่อนอื่นเรามาทำความรู้จัก solever แต่ะละตัวก่อน โดยเริ่มจากตัวแรก ode45 และ ode23 ที่ผมคัดรายละเอียดจากคู่มือมาสั้น ๆ ดังแสดงข้างล่าง

[t, y] = ode45 (fun, trange, init)

[t, y] = ode45 (fun, trange, init, ode_opt)

[t, y, te, ye, ie] = ode45 (…)

solution = ode45 (…)

ode45 (…)

Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Dormand-Prince method of order 4.

[t, y] = ode23 (fun, trange, init)

[t, y] = ode23 (fun, trange, init, ode_opt)

[t, y, te, ye, ie] = ode23 (…)

solution = ode23 (…)

ode23 (…)

Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Bogacki-Shampine method of order 3.

ode45 ใช้วิธี Dormand-Prince method อันดับ 4 ในขณะที่ ode23 ใช้วิธี Bogacki-Shampine อันดับ 3 เป็นตัวแก้ปัญหา ทั้งสองวิธีนี้เป็นวิธีที่พัฒนามาจาก Runge-Kutta อันดับ 4 ที่เป็นที่รู้จักกันทั่วไป รายละเอียดของวิธีการทั้งสองเคยเล่าไว้ใน Memoir ปีที่ ๑๐ ฉบับที่ ๑๔๔๘ วันพุธที่ ๔ ตุลาคม ๒๕๖๐ เรื่อง "การคำนวณเชิงตัวเลข (๑๘) เปรียบเทียบระเบียบวิธี Runge-Kutta" ความแตกต่างของสองวิธีการนี้ก็คือ ode23 จะให้คำตอบที่มีระดับความถูกต้องต่ำกว่า ode45 แต่การคำนวณทำได้รวดเร็วกว่า ตรงนี้ก็ขึ้นอยู่กับว่าผู้ใช้ต้องการคำตอบที่ถูกต้องขนาดไหน (เช่นคำตอบต้องถูกต้องถึงนัยสำคัญตำแหน่งที่ ๑๖ เลยหรือไม่) และความเร็วในการได้คำตอบเป็นเงื่อนไขที่สำคัญสำหรับการคำนวณหรือไม่ (เช่นในงานที่ต้องการนำเอาผลการคำนวณไปใช้ในการควบคุมกระบวนการ)

  

[t, y] = ode15s (fun, trange, y0)

[t, y] = ode15s (fun, trange, y0, ode_opt)

[t, y, te, ye, ie] = ode15s (…)

solution = ode15s (…)

ode15s (…)

Solve a set of stiff Ordinary Differential Equations (ODEs) or stiff semi-explicit index 1 Differential Algebraic Equations (DAEs).

ode15s uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.

solver ตัวที่สามคือ ode15s ที่มีรายละเอียดว่าใช้กับระบบสมการที่ "Stiff" หรือระบบสมการผสมพีชคณิต-อนุพันธ์ (Differential Algebraic Equations ที่ย่อว่า DAE) การแก้ปัญหาใช้วิธี Backward Differential Formula (ที่เป็นหนึ่งในกลุ่มวิธี multi-step) สำหรับผู้ที่สงสัยว่าสมการที่ "Stiff" นั้นเป็นอย่างไรสามารถอ่านได้ใน Memoir ปีที่ ๑๐ ฉบับที่ ๑๔๓๗ วันอังคารที่ ๑๒ กันยายน ๒๕๖๐ เรื่อง "การคำนวณเชิงตัวเลข (๑๗) เปรียบเทียบการแก้ Stiff equation ด้วยระเบียบวิธี Runge-Kutta และ Adam-Bashforth"

[t, y] = ode15i (fun, trange, y0, yp0)

[t, y] = ode15i (fun, trange, y0, yp0, ode_opt)

[t, y, te, ye, ie] = ode15i (…)

solution = ode15i (…)

ode15i (…)

Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or index 1 Differential Algebraic Equations (DAEs).

ode15i uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.

solver ตัวสุดท้ายคือ ode15i ที่ให้รายละเอียดว่าเป็นวิธี "fully-implicit" แต่มีการใช้วิธี Backward Differential Formula ร่วมด้วย

  

ทุกวิธีที่กล่าวมาข้างต้นใช้การคำนวณแบบ variable step คือระยะจากจุดตั้งต้นไปยังจุดถัดไปนั้นถูกกำหนดด้วยระดับความถูกต้องที่ต้องการของคำตอบ ดังนั้นผลการคำนวณที่ออกมาจะพบว่าระยะห่างระหว่างจุดข้อมูลนั้นไม่คงที่ แต่ความถูกต้องของคำตอบที่ได้นั้นจะคงที่ (เช่นสมมุติว่าต้องการคำตอบที่ถูกต้องในระดับนัยสำคัญ ๔ ตำแหน่ง ทุกคำตอบที่ได้มาก็จะมีคำตอบที่ถูกต้องที่ระดับนี้ โดยอาจจะมากได้บ้างแต่จะไม่น้อยกว่า) กล่าวคือในบริเวณที่คำตอบมีการเปลี่ยนแปลงอย่างรวดเร็ว อัลกอริทึมจะใช้ step size ที่แคบเพื่อรักษาความถูกต้องของผลการคำนวณ และในบริเวณที่คำตอบมีการเปลี่ยนแปลงช้า อัลกอริทึมก็จะใช้ step size ที่กว้างเพื่อทำให้การคำนวณดำเนินไปข้างหน้าได้อย่างรวดเร็ว

  

สมการที่นำมาใช้ทดสอบคือสมการอัตราการเกิดปฏิกิริยาของสาร A เมื่อเวลาใด ๆ ที่มีรูปแบบ dA/dt = -kA โดยที่เมื่อ t = 0, A = 1 และ k = 100 สมการนี้มีคำตอบคือ A = A0 exp (-kt) เมื่อ A₀ คือค่า A ที่เวลา t = 0

  

เพื่อทำการเปรียบเทียบผลการคำนวณด้วย solver ของ GNU octave จึงได้ทำการเขียนคำสั่งเพื่อหาค่า A ที่เวลา t ใด ๆ ด้วยเทคนิค implicit Euler โดยสมการที่ได้คือ A(i+1) = A(i)/(1+kΔt) โดยการคำนวณนี้ใช้ Δt = 0.001 คงที่ตลอดการคำนวณ

  

รูปที่ ๑ - ๓ คือชุดคำสั่งใช้ในการคำนวณ

  

รูปที่ ๑ ชุดคำสั่งสำหรับการคำนวณด้วย solver ode45, ode23 และ ode15s ทั้งสาม solver มีรูปแบบการเรียกใช้ฟังก์ชันที่เหมือนกัน เปลี่ยนเพียงแค่ชื่อของ solver ตัวที่ต้องการใช้ที่บรรทัดที่ 5 การจัดสมการจัดในรูปแบบ dy/dx = f(x,y) (ตรงบรรทัดที่ 4) ผลการคำนวณที่ได้ถูกบันทึกลงไฟล์ "result4.txt"

  

รูปที่ ๒ ชุดคำสั่งสำหรับการคำนวณด้วย solver ode15i ตัวนี้จะจัดสมการแตกต่างไปจากสามตัวแรก คือจัดในรูปแบบ dy/dx + f(x,y) = 0 และจำเป็นต้องมีการให้ค่า x และ dy/dx ที่เวลาเริ่มต้นด้วย ในกรณีของโจทย์ข้อที่ยกมาเป็นตัวอย่างนี้ค่า x ที่เวลาเริ่มต้นคือ A₀ = 1 และค่า dy/dx ที่เวลาเริ่มต้นคือ -kA₀ = -100 (ในที่นี้ใช้ค่า k₁ = 100) แต่ก็อาจใช้คำสั่ง decicช่วยในการหาค่าเริ่มต้นได้ด้วย สำหรับโจทย์ข้อนี้ทดลองดูแล้วไม่เห็นว่ามันจะต่างกัน ดังนั้นคำสั่งนี้สงสัยว่าน่าจะเหมาะสำหรับกรณีของระบบหลายสมการมากกว่า

คำตอบที่ได้จากการคำนวณแสดงไว้ในรูปที่ ๔ - ๘ (การคำนวณกระทำบนระบบปฏิบัติการ 64 บิต)

  

ผลการคำนวณด้วย solver ode45 ออกมาดีอย่างที่ควรเป็น เมื่อตรวจดูรายละเอียดค่าที่ solver คำนวณได้กับคำตอบที่ถูกต้องก็เห็นว่าใกล้เคียงกัน ความถี่ของจุดข้อมูลจะสูงในช่วงแรก (ช่วงที่ความเข้มข้นลดลงอย่างรวดเร็ว) เห็นได้จาก Δt มีขนาดเล็ก ก่อนที่ระยะห่างระหว่างจุดจะเพิ่มขึ้นในช่วงหลัง (ช่วงที่ความเข้มข้นเปลี่ยนแปลงอย่างช้า ๆ)

  

แต่ในกรณีของ ode 23 ให้ผลที่แตกต่างไปในช่วงท้าย (รูปที่ ๕) คือคำตอบในช่วงแรกออกมาดี แต่พอมาถึงบริเวณที่ความเข้มข้นเข้าใกล้ศูนย์มากกลับพบว่าค่ามีการแกว่งขึ้นสลับลงระหว่างบวกและลบ (ซึ่งไม่ถูกต้อง) สาเหตุหนึ่งตรงนี้เป็นผลมาจากการที่ค่า k ค่อนข้างจะสูง ทำให้สมการนี้มีความ stiff อยู่บ้าง

  

รูปที่ ๓ ชุดคำสั่งสำหรับการคำนวณด้วยเทคนิค Implicit Euler ที่ใช้สำหรับเปรียบเทียบคำตอบ

รูปที่ ๔ คำตอบที่ได้จากการคำนวณด้วย solver ode45 เมื่อใช้ k = 100

รูปที่ ๖ เป็นคำตอบที่ได้จากการใช้ solver ode15i จะเห็นได้ว่าในกรณีของตัวอย่างนี้คำตอบที่ได้อยู่ในระดับเดียวกันกับเมื่อใช้ ode45 เพียงแต่การใช้ ode15i จะใช้ระยะ Δt ที่แคบกว่า และเมื่อเปลี่ยนไปใช้ ode15i กลับพบพฤติกรรมที่ไม่คาดว่าจะพบคือ คำตอบที่ได้มีการแกว่งขึ้นลงในช่วงท้าย (เห็นได้ชัดจากการที่เครื่องหมายของคำตอบมีการสลับบวกลบ) ซึ่งไม่น่าจะเป็นพฤติกรรมที่เกิดขึ้นกับเทคนิคตระกูล implicit ดังผลที่แสดงในรูปที่ ๙ กราฟในรูปที่ ๙ ที่เห็นไม่แนบชิดกับคำตอบที่ถูกต้องเป็นเพราะ Δt ที่ใช้นั้นกว้างเกินไป ถ้าลด Δt ให้มีขนาดเล็กลงก็จะได้ผลที่ใกล้กันเข้าไปอีก

  

รูปที่ ๕ คำตอบที่ได้จากการคำนวณด้วย solver ode23 เมื่อใช้ k = 100

  

รูปที่ ๖ คำตอบที่ได้จากการคำนวณด้วย solver ode15s เมื่อใช้ k = 100

การที่ ode15i ที่ซอร์ฟแวร์อ้างว่าใช้เทคนิค fully implicit ในการแก้ปัญหาแล้วคำตอบมีการแกว่งนี้เป็นเรื่องแปลก ทั้ง ๆ ที่เทคนิคตระกูลนี้มีความสามารถในการควบคุมการสะสมของ error ที่ดี เรื่องนี้คงต้องค่อย ๆ ศึกษากันต่อไป

  

รูปที่ ๗ คำตอบที่ได้จากการคำนวณด้วย solver ode15i เมื่อใช้ k = 100

  

รูปที่ ๘ คำตอบที่ได้จากการคำนวณด้วยเทคนิค Implicit Euler เมื่อใช้ k = 100

ในตอนหน้าเราจะมาดูกันว่า ถ้าเพิ่มค่า ให้สูงขึ้นไปอีก k ผลการคำนวณจะออกมาเป็นอย่างไร