การแก้สมการ Differential ด้วย MATLAB

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

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

ตัวอย่างเช่น ผมมีสมการ differential ที่ต้องการแก้อยู่ 3 สมการ ดังนี้

dy1/dt = y2*y3

dy2/dt = -y1*y3
dy3/dt = -0.51*y1*y2

โดยมีค่าเริ่มต้น (initial condition) ดังนี้

y1(0) = 0

y2(0) = 1

y3(0) = 1

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

ง่ายมากเลยครับ เราก็แค่ลองเปลี่ยนค่าเริ่มต้นดูหลายๆ ค่า ถ้าเป็นค่าเริ่มต้นที่ถูกต้อง คำตอบจะออกมาคล้ายๆ กัน (แต่ไม่เหมือนกันเป๊ะนะครับ แค่คล้ายๆ กัน) แต่ถ้าค่าเริ่มต้นนั้นใช้งานไม่ได้ คำตอบหรือกราฟที่ได้ จะต่างกันโดยสิ้นเชิงครับ

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

[ X , Y ] = ode45(@eqfcn, x_span , y_init)

โดยที่

X และ Y คือคำตอบ

eqfcn คือ ฟังก์ชันของสมการ (เราจะต้องเขียนสมการ diff แยกไว้ในฟังก์ชันต่างหาก)

x_span คือ ช่วงค่าของ x ที่เราอยากดู

y_init คือ เริ่มเริ่มต้นของ y

เอาละ เรามาดูตัวอย่างกันเลยครับ

อันดับแรก สร้างไฟล์ที่ชื่อว่า 'simple1.m' ขึ้นมาก่อนครับ ซึ่งในนี้เราจะสร้างฟังก์ชันชื่อว่า 'simple1' เรื่องนี้สำคัญนะครับ และหลายๆ คนมักจะทำผิดตรงจุดนี้ ทุกคนต้องจำไว้ให้ดีนะครับ เพราะนี้คือกฎการเขียนฟังก์ชันข้อที่ 1 ชื่อไฟล์ และชื่อฟังก์ชัน ต้องตรงกันเสมอ (ยกเว้น private function)

function dy = simple1(t,y)
dy = zeros(3,1);
y1 = y(1);
y2 = y(2);
y3 = y(3);
dy(1) = y2*y3;
dy(2) = -y1*y3;
dy(3) = -0.51*y1*y2;
end

สังเกตุดูนะครับ ตัวแปร t ซึ่งเป็น input ของฟังก์ชัน ไม่มีการนำไปใช้งานในฟังก์ชันเลย แต่เราต้องใส่เอาไว้ครับ เพราะมันคือส่วนของ dt

ส่วนตัวแปร y1, y2, y3 ผมเขียนขึ้นมาเพื่อให้เห็นได้ง่ายๆ เท่านั้นครับ สำหรับคนที่ยังไม่เก่งเรื่อง อาเรย์ เพราะเราจะเห็นว่าตัวแปร y ซึ่งเป็น input มีแค่ตัวเดียวเท่านั้น แต่สมการของเรามีถึง 3 สมการ หลายๆ คนอาจจะสงสัยว่า แล้วมันจะใส่ค่าทั้ง 3 สมการเข้ามาได้ไง

คำตอบก็คือ MATLAB จะใส่ค่าเข้ามาในรูปแบบของอาเรย์ครับ โดยมีจำนวณเท่ากับค่า initial ที่เรากำหนด ดังนั้นเราก็ต้องเข้าใจว่า y(1) ก็คือ y1 นะ ส่วน y(2) ก็คือ y2 และ y(3) ก็คือ y3 ตามลำดับ แต่ถ้าใครที่เข้าใจเรื่องอาเรย์อยู่แล้ว จะแทน y1 ในสมการด้วย y(1) เลยก็ได้ และก็ไม่จำเป็นต้องเขียน 3 บรรทัดนี้ครับ

y1 = y(1);
y2 = y(2);
y3 = y(3);

โปรแกรมเราก็จะดูสั้นลงไปเยอะ

จากนั้นเรามาลองมาแก้สมการพวกนี้ด้วยคำสั่ง ode45 ดูครับ

t_span = [0 12];      %Time period
Y_init = [0 1 1];     %Initial condition
[T1,Y1] = ode45(@simple1,t_span,Y_init);
plot(T1,Y1(:,1),'-', T1,Y1(:,2),'-.', T1,Y1(:,3),':');

คำสั่งชุดนี้ ให้เขียนแยกเอาไว้ในไฟล์อื่นนะครับ หรือเขียนใน command window ก็ได้ แต่อย่าเขียนรวมกันไว้ในไฟล์ simple1 (แต่ถ้าผู้หัด เขียน private function เป็น จะเขียนไว้ด้วยกันก็ได้)

ในโค้ด ผมกำหนดค่า t_span จาก 0 ถึง 12 หมายความว่าเราจะดูค่าในช่วงเวลา t = 0 ถึง 12 นะครับ
ส่วน y_init ก็คือค่าเริ่มต้นที่ผมได้อธิบาย ซึ่งมีทั้งหมด 3 ค่า จากนั้นกดปุ่ม RUN ได้เลยครับ

ผลที่ได้

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