บทที่ 8 การแก้สมการอนุพันธ์

บทที่  8

การแก้สมการอนุพันธ์

8.1 ส่วนนำทาง

            เอกสารนี้ออกแบบเพื่อให้ผู้ใช้ได้คุ้นเคยกับการแก้สมการอนุพันธ์สามัญ(ordinary differential equations:ode) โดยใช้ขั้นตอนวีธี  ode ไซแลบ

8.2 สมการอนุพันธ์

            8.2.1 สเกลล่า ODE’s

            สมมุติว่าต้องการหาคำตอบ u(t) ของสมการอนุพันธ์สามัญดังนี้

                                                                                (8.1)

            สมการนี้เป็นรูปแบบอย่างง่ายสำหรับการแพร่กระจายเชื้อโรค  พิจารณาฝูงสัตว์ P และกำหนดให้ u(t) เป็นสัดส่วนของสัตว์ที่ติดเชื้อโรคหลังจากเวลาผ่านไป t วัน ดังนั้น 1 − u(t) คือสัดส่วนของสัตว์ที่ไม่ติดเชื้อ  การติดเชื้อใหม่จะเกิดขึ้นเมื่อสัตว์ที่ติดเชื้อกับสัตว์ที่ไม่ติดเชื้อมาพบปะกัน จำนวนครั้งของการพบปะกันเป็นสัดส่วนตรงกับ  u(t) (1 − u(t)).  ดังนั้นแบบรูป(model) สำหรับการแพร่ขยายเชื้อโรค คือจำนวนที่ติดเชื้อใหม่ต่อวัน  du/dt, เป็นสัดส่วนกับ u(t) (1 − u(t)). ตามที่ได้ในสมการที่ 8.1  ซึ่งตัวคงที่ k ขึ้นอยู่กับความหนาแน่นของสัตว์ และการติดเชื้อ  สำหรับความแน่ชัด สมมุติว่า P เป็นสัตว์ 10000 ตัวเขียนได้ว่า P= 10000  และ k = .2. เงื่อนไขขั้นต้นที่เวลา t = 0 มีสัตว์ที่ติดเชื้อเพียงตัวเดียว เพิ่มเงื่อนไขเริ่มแรกได้ดังนี้

            u(0) = 1/10000.                                                                                      (8.2)

            เมื่อต้องการหาว่าสัตว์แต่ละตัวจะติดเชื้อเท่าใดในอีก 100 วันข้างหน้า ทำได้โดยการพล็อต u(t)  จาก t = 0 ถึง t = 100. เพื่อแก้ปัญหานี้จำต้องกำหนดฟังก์ชันเพื่อใช้ในการคำนวณทางด้านซ้ายของสมการ 8.1  จากนั้นจึงใช้โปรแกรม ode ของไซแลบ

การสร้างฟังก์ชันไซแลบซึ่งระบุทางด้านขวาของสมการอนุพันธ์  ให้แน่ใจว่าได้เรียกลำดับอาร์กิวเมนต์ที่ถูกต้เอง  ดังเช่น

function udot=f(t,u)

    udot = k * u * ( 1 - u );

endfunction

(แม้ว่าไม่ได้ใช้ในการคำนวณ udot, ฟังก์ชัน odeไซแลบ  ฟังก์ชัน ode ต้องการฟังก์ชัน f เพื่อส่งผ่านมายัง ode จะต้องมี 2 อินพุตอาร์กิวเมนต์  (ในลำดับที่ถูกต้อง)  ให้ขอความช่วยโดยใช้ help ode สำหรับข้อสนเทศเพิ่มเติม
            กำหนดปัญหาขึ้นสมการอนุพันธ์สามารถแก้ได้เพียงบรรทัดเดียวในไซแลบโดยใช้ ode

t=0:100;

k = 0.2;

u = ode( 1/10000, 0.0, t, f );

(xbasc()); plot2d( t,u) // ถ้าใส่วงเล็บข้างหน้าแล้ว error

หมายเหตุ ค่าของ k หาได้สำหรับ ode และ f  ฟังก์ชัน ode มีพารามีเตอร์ที่สามารถปรับเพื่อควบคุมความละเอียดในการคำนวณ  ส่วนมากแล้วไม่จำเป็นสำหรับเพียงการพล็อตกราฟคำตอบของปัญหาธรรมดาสามัญ   เป็นแนวคิดที่ดีเสมอที่จะตรวจสอบคำตอบเดิมโดยแก้ปัญหาจากความต้องการที่ให้ค่าความละเอียดถูกต้องมากขึ้น  ตัวอย่างเช่น สามารถที่จะใช้ ode   ดังนั้นข้อผิดพลาดที่ประมาณอย่างสัมบูรณ์ ที่สัมพันธ์ที่จะจำกัดให้น้อยกว่า 1 × 10−9.  ทำเช่นนี้ได้ด้วยคำสั่ง

uu = ode( 1/10000, 0.0, t, rtol=1e-9, atol=1e-9, f );

plot2d(t,uu,style=-1)

จะเห็นได้ทันที่ว่าผลเฉลยที่คำนวณได้ในการเรียกใช้ครั้งแรกกับ ode วางอยู่บนสุดได้เคิร์ฟที่ละเอียดแม่นยำกว่า  นี่เป็น 2 ผลเฉลยที่ให้ความละเอียดเชิงกราฟิกส์เหมือนกัน  สำหรับตัวอย่างนี้รูปแบบหรือโมเดลเป็นไปได้ที่จะหาสูตรวิเคราะห์สำหรับu(t).  จริงแล้ว

u(t) = 1/(1 + c exp(−kt))   ที่ซึ่ง
            c = P − 1.                                                                                              (8.3)

            โดยเพิ่มสูตรนี้เข้าไปในกราฟข้างบนเป็นการตรวจสอบสุดท้าย

t = 0:100 ;

c = 10000-1 ; k = .2 ;

function u=exactu(t)

     u = 1 ./ ( 1 + c*exp(-k*t));

endfunction

xbasc(); plot2d(t,u); fplot2d(t,exactu,style=-1 )

            สูตรที่ชัดเจนให้เส้นเคิร์ฟเช่นเดียวกับผลเฉลยเชิงตัวเลข จริงแล้วข้อแตกต่างระหว่างผลเฉลยทั้งสองมากที่สุดที่ 2.5×10−7. ดังนั้นจึงมั่นใจได้ว่าการใช้ ode ในการแก้ปัญหาข้างล่างต่อไป

แบบฝึกหัดที่13.

            1. จงตรวจสอบว่า (8.3) คือผลเฉลยของ (8.1).  นั่นคือ เมื่อทำให้เล็กย่อยลง(differentiate) สูตรนี้สำหรับ u(t) จะได้ทางขวาของ (8.1).  และสูตรนี้สอดคล้องกับ (8.2)

2. จากกราฟ u(t), ประมาณการได้ว่าเป็นจำนวนกี่วันที่ต้องการที่จำทำให้ผลเมืองติดเชื้อได้ครึ่งหนึ่ง

แนะ:สามารถตรวจสอบการพล็อตเพื่อให้ได้ผลเฉลยจำนวนวันที่ใกล้เคียงมากที่สุด  คำสั่ง xgrid จะทำให้ง่ายในการตอบคำถาม  จึงควรที่จะสร้างฟังก์ชันหนึ่งโดยใช้ผลจากขั้นตอนกระบวนการ ซึ่งคืนผลเฉลย ode เป็นลบ0.5 ที่เวลาเฉพาะ t ใดๆ  แล้วใช้ fsolve เพื่อแก้ปัญหา

            3. สมมุติว่าตอนนี้หลังจากการติดเชื้อบันทึกไว้ที่วันที่  t = 0,  โปรแกรมการฉีดวัคซีนเริ่มจากวันที่ 5 และ 200  ของสัตว์ที่ยังไม่ติดเชื้อ ได้รับการฉีดวัคซีนต่อวัน  หลังจาก t วัน ประชากรของสัตว์ที่ฉีดวัคซีนต่อวันคือ  200(t−5)+/10000.  จากนี้ (t−5)+ means max(t−5, 0). สัตว์ที่ได้รับวัคซีนเหล่านี้ไม่ติดเชื้อ ส่วนที่ติดเชื้อใหม่เกิดขึ้นเมื่อสัตว์ที่ไม่ติดเชื้อ สัตว์ที่ไม่ได้รับวัคซีนได้พบร่วมกับสัตว์ที่ติดเชื้อ  ดังนั้นอัตราการเพิ่มของ u ตอนนี้สร้างแบบรูปโดยสมการอนุพันธ์ดังนี้

du

dt      = k u(t) (1 − u(t) − .02 (t − 5)+)+ ;

            จากนี้ให้หาว่ามีสัตว์เท่าใดที่ไม่ติดเชื้อหลังจาก 100 วันไปแล้วที่ได้รับโปรแกรมวัคซีน

            4. มีสัตว์จำนวนเท่าใดที่ไม่ติดเชื้อ ถ้าเฉพาะสัตว์ 100 ตัวฉีดวัคซีนในแต่ละวันจากวันที่ 5  ,มีสัตว์เท่าใดที่จะปลอดภัยถ้า สัตว์ 100 ตัวต่อวันที่ได้รับการฉีวัคซีนากวันที่ 0?  ปัญหาแบบรูปอย่างง่ายนี้ มีสูตรที่ใช้หาผลเฉลยที่แน่ชัด  อย่างไรก็ตามในไซแลบนั้นจะง่ายมากที่จะให้ได้ผลเฉลยเป็นจำนวนตัวเลข และทำให้เสร็จสิ้นลงไม้ว่าจะเป็นปัญหาที่ซับซ้อนก็ตาม  เมื่อไม่มีสูตรที่ง่ายอื่นใดสำหรับหาผลเฉลย  การใช้ไซแลบสามารถสำรวจแบบรูปได้อย่างง่ายดายโดยการเปลี่ยนค่าพารามีเตอร์ และการพล็อตของผลเฉลย

8.2.2 ODE ลำดับสอง

            ตามตัวอย่างที่สองนี้ พิจารณาสมการอนุพันธ์ธรรมดาทั่วไป

                                                                            (8.4)

            u(0) = π/4                                                                          (8.5)

นี่เป็นตัวอย่างที่ซับซ้อนมากขึ้น เพราะว่าปรากฏของอนุพันธ์ของ u ลำดับ2  ไม่เป็นเพียงอนุพันธ์บำดับ 1 . (เรียกว่าว่าเป็น ODE ลำดับ2)  ส่วนแก้สมการ ODE ไซแลบไม่สามารถที่จะใช้อนุพันธ์ลำดับสองได้โดยตรง เหมือนกับที่ทำงานได้กับอนุพันธ์ลำดับ1  อย่างไรก็ตามก็ยังมีกลวิธีที่สามารถนำมาใช้  ถ้าหากติดติดตามทั้ง  u(t) และ du/dt, สำหรับอนุพันธ์ลำดับ2  ของ u(t) ก็คิดให้เป็นอนุพันธ์ลำดับแรกของ du/dt.  ที่เป็นแบบฉบับมากขึ้น กำหนด  v(t) เป็นค่าของ du/dt.  แล้วสมการที่ (8.2) เทียบเท่ากับ 2 สมการดังนี้

และ

แม้แต่ให้รูปสะสวยขึ้น โดยกำหนด  x(t) เป็นเว็คเตอร์ ที่มี 2 องค์ประกอบ

        

และดังนั้น

           

ตอนนี้อยู่ในรูป

ดังนั้นการแก้สมการ (8.2) และ  plot u(t)  กับ t   โดยเริ่มโดยนิยามฟังก์ชัน F

function xdot = pendulum(t,x)

xdot = [ x(2) ; -sin(x(1))];

endfunction

Then we use the commands

t = 0:.25:20;

x=ode([%pi/4;0], 0.0, t, pendulum );

xbasc(); plot2d(t, x(1,:))

            เพื่อที่จะพล็อตองค์ประกอบแรกของผลเฉลย  สังเกตได้ว่า ode เก็บองค์ประกอบที่1 และ 2 ของผลเฉลยในแถวที่1 และ 2 ของเอ้าพุต  แน่นอนว่าสามารถพล็อตองค์ประกอบที่ 2 ของผลเฉลยและ ผิวหน้า(phase)ของระนาบเส้นทางการเคลื่อนที่ด้วยคำสั่งดังนี้

xbasc(); plot2d(t, x(2,:))

xinit(); xbasc(); plot2d(x(1,:), x(2,:))

            คำสั่ง xinit() ทำอะไร?

            ความเร็วสนามในด้านของระนาบสามารถที่จะวาดโดยใช้คำสั่ง fchamp ดังนี้

fchamp(pendulum,0.0,-1:0.1:1,-1:0.1:1)

Exercise 14.  ในสมการ(8.2) , u(t)  คือมุมที่แก่วงของลูกตุ้มนาฬิกาที่ทำกับแกนในแนวตั้งที่เวลา t ถ้าลูกตุ้มนาฬิกามีความยาว l น้ำหนักของลูกตุ้มนาฬิกาที่ตำแหน่ง  (x, y) = (sin(u(t)),−cos(u(t)). (ลูกตุ้มนาฬิกาแกว่งจากจุดกำเนิด (origin) ที่ตำแหน่ง(0,−1) เมื่ออยู่ในสภาพนิ่ง)

1. จากกราฟ ให้หาคาบของลูกตุ้มนาฬิกา

2. บ่อยครั้งเมื่อลูกตุ้มแกว่งเป็นมุมเล็กๆ ตามสมการ (8.2) ประมาณการโดย

           

u(0) = π/4 ,

            ผลเฉลยกับการประมาณเชิงเส้น และผลเฉลยที่แน่ชัด แตกต่างกันอย่างไร  คาบเวลาเท่ากันหรือไม่3. ลูกตุ้มนาฬิกาแบบบังคับสามารถจัดแบบรูปโดยสมการ

                                                                           (8.6)

ก่อให้ได้คำตอบของสมการนี้สำหรับr ω = 1/2 และ A = 0.2  และ A = 1.5 ด้วยเงื่อนไขเริ่มแรก

u(0) = 100, u′(0) = 0  และ u(0) = 101, u′(0) = 0.

ความคิดเห็นเกี่ยวกับความไวของคำตอบของ ode สำหรับค่าที่แตกต่างของ A .