#!/usr/bin/env python # coding: utf-8 # # พาราโบลาตัดแกน $x$ ที่ไหน # # แนะนำให้ดู [Visual introduction to parabolas](https://www.youtube.com/watch?v=BGz3pkoGPag) # # พาราโบลาทุกอันมีรูปทรงเดียวกัน ต่างกันที่เราขยับ หมุน และซูมเข้าออกอย่างไร (เหมือนกับวงกลมทุกอันคล้ายกันหมด แค่ซูมเข้าออก) แนะนำให้ดู [There is only One True Parabola](https://www.youtube.com/watch?v=hoh4TmPzu1w) # # ปกติเราจะเขียนสมการพาราโบลาในรูป $y(x) = a x^2 + b x + c$ โดยที่ $a$ ไม่เป็น 0 # # ถ้าเราจะหาว่าพาราโบลาตัดแกน $x$ ที่ไหน เราก็ต้องหาค่า $x$ ที่ทำให้ $y(x) = 0$ ก็คือการแก้สมการ $a x^2 + b x + c = 0$ นั่นเอง # # เราแก้สมการได้ด้วยวิธีที่เรียกว่า [completing the square:](https://en.wikipedia.org/wiki/Completing_the_square) คือจัดรูปให้เป็นกำลังสองของอะไรบางอย่างเท่ากับอะไรบางอย่างอีกก้อนหนึ่ง: # # $a x^2 + b x + c = 0$ # # $ x^2 + \frac{b}{a} x + \frac{c}{a} = 0$ # # $ x^2 + 2 (\frac{b}{2 a}) x + \frac{c}{a} = 0$ # # $ x^2 + 2 (\frac{b}{2 a}) x + (\frac{b}{2 a})^2 -(\frac{b}{2 a})^2 + \frac{c}{a} = 0$ # # $ (x + (\frac{b}{2 a}))^2 -(\frac{b}{2 a})^2 + \frac{c}{a} = 0$ # # $ (x + (\frac{b}{2 a}))^2 = (\frac{b}{2 a})^2 - \frac{c}{a} $ # # $ (x + (\frac{b}{2 a}))^2 = \frac{b^2 - 4 a c}{4 a^2}$ # # $ x + (\frac{b}{2 a}) = \pm \sqrt{\frac{b^2 - 4 a c}{4 a^2}}$ # # $ x = -(\frac{b}{2 a}) \pm \sqrt{\frac{b^2 - 4 a c}{4 a^2}}$ # # $ x = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}$ # # ดังนั้น ค่า $ x = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}$ จะเป็นค่า $x$ ที่ทำให้พาราโบลา $a x^2 + b x + c$ ตัดแกน $x$ # # # # In[1]: # เราสร้างฟังก์ชั่นเพื่อคำนวณตามความสัมพันธ์ข้างบน import math def zeroes_of_parabola(a,b,c): """หาค่า x ที่ทำให้ a x**2 + b x + c == 0 ค่า x ที่หาได้คือจุดที่พาราโบลาตัดแกน x a ต้องไม่เป็น 0 """ if a == 0: # ถ้า a เป็น 0 จะไม่่คำนวณอะไรให้ return None if b**2 - 4*a*c < 0: #ถ้าสิ่งที่อยู่ในสแควร์รูทเป็นลบก็จะไม่คำนวณอะไรให้ return None if b**2 - 4*a*c == 0: #ถ้า b**2 เท่ากับ 4*a*c พอดีจะมีคำตอบเดียว return [-b/(2*a)] return [(-b + math.sqrt(b**2-4*a*c))/(2*a),\ (-b - math.sqrt(b**2-4*a*c))/(2*a)] #ไม่งั้นจะมีสองคำตอบ # In[2]: # -x^2 + 4 = 0 มีคำตอบที่ x = ± 2 zeroes_of_parabola(-1,0,4) # In[3]: # x^2 - 2x + 1 มีคำตอบที่ x = 1 zeroes_of_parabola(1,-2,1) # In[4]: # x^2 - 2x + 4 ไม่มีคำตอบที่เป็นจำนวนจริงเพราะ b^2 - 4ac น้อยกว่า 0 zeroes_of_parabola(1,-2,4) # # ถ้าเราไม่มีสูตรหาคำตอบ เราอาจใช้วิธีประมาณเช่น Bisection Method # # สมการของพาราโบลาที่มี $x$ ยกกำลังสูงสุดเท่ากับ $x^2$ เรียกว่าสมการควอดราติก ([quadratic equation](https://en.wikipedia.org/wiki/Quadratic_equation)) # # เรามีสูตร $ x = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}$ สำหรับแก้สมการควอดราติก $a x^2 + b x + c = 0$ # # ถ้าสมการมี $x^3$ เป็นกำลังสูงสุดเราเรียกว่าคิวบิก ([cubic equation](https://en.wikipedia.org/wiki/Cubic_function)) ถ้ามี $x^4$ เป็นกำลังสูงสุดเราเรียกว่าควอติก ([quartic equation](https://en.wikipedia.org/wiki/Quartic_function)) สมการกำลังห้าเรียกว่าควินติก ([quintic equation](https://en.wikipedia.org/wiki/Quintic_function)) # # เมื่อประมาณสี่ห้าร้อยปีมาแล้วนักคณิตศาสตร์พบสูตรแก้สมการคิวบิกและควอติก กล่าวคือมีคนค้นพบสูตรสำหรับแก้สมการ $a x^3 + b x^2 + c x + d = 0$ และ $a x^4 + b x^3 + c x^2 + d x + e = 0$ แม้ว่าสูตรจะยุ่งยากและยาวกว่าสูตรของควอดราติกมากก็ตาม เวลาผ่านไปอีกเกือบๆ 300 ปีก่อนที่จะมี[คนค้นพบ](https://en.wikipedia.org/wiki/Abel–Ruffini_theorem)ว่าสมการที่มีกำลังตั้งแต่ 5 ขึ้นไปจะไม่มีสูตรทั่วไปที่จะแก้สมการเหล่านั้น # # สำหรับสมการที่ไม่มีสูตรสำหรับแก้ เราสามารถหาคำตอบโดยการประมาณได้หลายวิธี วิธีแรกที่จะใช้วันนี้เรียกว่าวิธีแบ่งครึ่งหาคำตอบ ([bisection method](https://en.wikipedia.org/wiki/Bisection_method)) โดยเราจะอาศัยเทคนิควาดกราฟเพื่อดูก่อนว่าคำตอบของสมการ $f(x) = 0$ อยู่แถวๆไหน แล้วกำหนดช่วง $(xmin, xmax)$ ที่คำตอบเราอยู่แถวๆนั้นโดยระวังว่าเครื่องหมายของ $f(xmin)$ และ $f(xmax)$ ต้องตรงกันข้ามกัน คือถ้าฟังก์ชั่น $f$ มีการเปลี่ยนเครื่องหมายจากบวกไปลบหรือลบไปบวกในช่วง $(xmin, xmax)$ ก็แสดงว่าค่าฟังก์ชั่นต้องเป็น 0 ที่ไหนสักแห่งในช่วง $(xmin, xmax)$ นั้น จากนั้นเราก็จะแบ่งครึ่งช่วง $(xmin, xmax)$ เป็นสองช่วงเล็กๆ แล้วดูว่าฟังก์ชั่น $f$ เปลี่ยนเครื่องหมายในช่วงไหน ถ้าเปลี่ยนในช่วงไหนก็แสดงว่าคำตอบ $f(x) = 0$ อยู่ในช่วงนั้น เราสามารถแบ่งครึ่งช่วงอย่างนี้ไปเรื่อยๆให้ช่วงเล็กลงมากๆได้ พอช่วงมีขนาดเล็กพอเราก็บอกว่าคำตอบ $f(x) = 0$ อยู่ในช่วงเล็กๆนั้น อาจจะใช้ค่าตรงกลางช่วง (พอใช้ได้) หรือเทียบบัญญัติไตรยางค์ระหว่างค่าฟังก์ชั่นที่ขอบทั้งสองของช่วงก็ได้ (ได้คำตอบดีกว่า) # In[5]: # ตัวอย่างฟังก์ชั่นกำลังห้าของเรา def quintic(x): return x**5 + x**4 - x**3 + 2*x**2 - 5*x + 6 # กราฟฟังก์ชั่นกำลังห้าของเรา ($x^5+x^4-x^3+2x^2 - 5x + 6)$ จะเห็นได้ว่าค่าฟังก์ชั่นเป็นศูนย์ (ตัดแกน $x$) ในช่วง (-3, -2): # # ![Screen%20Shot%202562-06-29%20at%2011.09.02.png](attachment:Screen%20Shot%202562-06-29%20at%2011.09.02.png) # In[6]: # ฟังก์ชั่น bisection method เวอร์ชั่นแรก # พิมพ์ว่า xmin, xmax มีขนาดเท่าไร ผลคูณ f(xmin)*f(xmax) เป็นเท่าไร # หยุดทำงานเมื่อขนาด xmin ห่างจาก xmax ไม่เกิน 0.0001 # bisection มีการเรียกตัวเองด้วย (recursion) def bisection1(f, xmin, xmax): "พยายามหาค่า x ที่ทำให้ f(x)==0, xmin และ xmax คือช่วงที่เราเดาว่าคำตอบอยู่ในนั้น" if xmin > xmax: # จัดการให้ xmin น้อยกว่า xmax เสมอ xmin, xmax = xmax, xmin print(f"xmin, xmax = {xmin}, {xmax}\t f(xmin)*f(xmax) = {f(xmin)*f(xmax):.5f}") if xmax-xmin < 0.0001: print(f'Found a solution: {(xmax+xmin)/2}') return( (xmin+xmax)/2 ) if f(xmin) == 0: #ถ้า xmin หรือ xmax ทำให้ค่าฟังก์ชั่นเป็นศูนย์ก็ตอบ xmin หรือ xmax เลย return xmin if f(xmax) == 0: return xmax if f(xmin)*f(xmax) > 0: #ถ้าฟังก์ชั่น f ไม่เปลี่ยนเครื่องหมายจากบวกไปลบหรือลบไปบวก #ระหว่าง xmin และ xmax ก็ควรไปเดาใหม่ว่า xmin, xmax คืออะไร print("Bad xmin, xmax") return None xmid = (xmin + xmax)/2 #หาจุดกลางระหว่าง xmin และ xmax if f(xmin)*f(xmid) < 0: #ถ้าฟังก์ชั่นเปลี่ยนเครื่องหมายระหว่าง xmin กับจุดกลาง ก็ไปหาต่อในช่วงนี้ x = bisection1(f,xmin, xmid) else: #ไม่งั้นก็หาต่อในช่วงจุดกลางถึง xmax x = bisection1(f,xmid, xmax) return x # In[7]: # พยายามหาค่า x ที่ทำให้ฟังก์ชั่นข้างบนของเรามีค่าเท่ากับ 0 โดยเริ่มเดาว่า x อยู่ระหว่าง -3 และ -2 x = bisection1(quintic,-3,-2) # In[8]: #ตรวจสอบว่าค่าที่หามาได้ทำให้ฟังก์ชั่นใกล้ๆศูนย์ไหม print(f"x = {x:.3f}, f(x) = {quintic(x):.3f}") # In[9]: # เวอร์ชั่นที่ 2 กำหนดว่าจะหยุดหาเมื่อ xmin ห่างจาก xmax ไม่เกิน tolerance ที่เราตั้งค่าไว้ # ฟังก์ชั่นเรียกตัวเองด้วย (recursion) def bisection2(f, xmin, xmax): "พยายามหาค่า x ที่ทำให้ f(x)==0, xmin และ xmax คือช่วงที่เดาว่าคำตอบอยู่ในนั้น" tolerance = 1e-6 #ตั้งค่่าไว้ว่าให้หยุดหาคำตอบเมื่อ xmin และ xmax ห่างกันไม่เกิน tolerance if xmin > xmax: # จัดการให้ xmin น้อยกว่า xmax เสมอ xmin, xmax = xmax, xmin if f(xmin)*f(xmax) > 0: print("Bad xmin, xmax") #ถ้าฟังก์ชั่น f ไม่เปลี่ยนเครื่องหมายจากบวกไปลบหรือลบไปบวก #ระหว่าง xmin และ xmax ก็ควรไปเดาใหม่ว่า xmin, xmax คืออะไร return None if xmax-xmin < tolerance: return( (xmin+xmax)/2 ) if f(xmin) == 0: #ถ้า xmin หรือ xmax ทำให้ค่าฟังก์ชั่นเป็นศูนย์ก็ตอบ xmin หรือ xmax เลย return xmin if f(xmax) == 0: return xmax xmid = (xmin + xmax)/2 #หาจุดกลางระหว่าง xmin และ xmax if f(xmin)*f(xmid) < 0: #ถ้าฟังก์ชั่นเปลี่ยนเครื่องหมายระหว่าง xmin กับจุดกลาง ก็ไปหาต่อในช่วงนี้ x = bisection2(f,xmin, xmid) else: #ไม่งั้นก็หาต่อในช่วงจุดกลางถึง xmax x = bisection2(f,xmid, xmax) return x # In[10]: x = bisection2(quintic,-3,-2) print(f"Value at x = {x:.6f} is {quintic(x):.6f}") # In[11]: # เวอร์ชั่น 3 ฟังก์ชั่นไม่เรียกตัวเองแล้ว (non-recursive) # ยังพิมพ์ดูว่า xmin, xmax, f(xmin)*f(xmax) เป็นเท่าไร # คำตอบที่ได้จะเป็นตรงกลางของ xmin และ xmax def bisection3(f,xmin, xmax): "พยายามหาค่า x ที่ทำให้ f(x)==0, xmin และ xmax คือช่วงที่เดาว่าคำตอบอยู่ในนั้น" tolerance = 1e-6 #ตั้งค่า tolerance ไว้ให้หยุดทำงาน #จะหยุดทำงานเมื่อ xmin ห่างจาก xmax น้อยกว่า tolerance if xmin > xmax: # จัดการให้ xmin น้อยกว่า xmax เสมอ xmin, xmax = xmax, xmin if f(xmin)*f(xmax) > 0: print("Bad xmin, xmax") #ถ้าฟังก์ชั่น f ไม่เปลี่ยนเครื่องหมายจากบวกไปลบหรือลบไปบวก #ระหว่าง xmin และ xmax ก็ควรไปเดาใหม่ว่า xmin, xmax คืออะไร return None while xmax-xmin >= tolerance: #จะหยุดทำงานเมื่อ xmin ห่างจาก xmax น้อยกว่า tolerance print(f"xmin, xmax = {xmin}, {xmax}\t f(xmin)*f(xmax) = {f(xmin)*f(xmax):.5f}") if f(xmin) == 0: return xmin if f(xmax) == 0: return xmax xmid = (xmin + xmax)/2 #หาจุดกลางระหว่าง xmin และ xmax if f(xmin)*f(xmid) < 0: #ถ้าฟังก์ชั่นเปลี่ยนเครื่องหมายระหว่าง xmin กับจุดกลาง ก็ไปหาต่อในช่วงนี้ xmin, xmax = xmin, xmid else: #ไม่งั้นก็หาต่อในช่วงจุดกลางถึง xmax xmin, xmax = xmid, xmax return (xmin+xmax)/2 #คำตอบที่ได้จะเป็นตรงกลางของ xmin และ xmax ที่ห่างกันน้อยกว่า tolerance # In[12]: x = bisection3(quintic,-3,-2) print(f"Value at x = {x} is {quintic(x)}") # In[13]: # เวอร์ชั่น 4 ฟังก์ชั่นไม่เรียกตัวเองแล้ว (non-recursive) # ไม่พิมพ์ดูว่า xmin, xmax, f(xmin)*f(xmax) เป็นเท่าไรแล้ว # คำตอบที่ได้จะเทียบบัญญัติไตรยางค์ระหว่างจุด (xmin, f(xmin)) และ จุด (xmax, f(xmax)) # ซื่งน่าจะได้คำตอบที่ใกล้ความจริงมากกว่าแบบตอบตรงกลางระหว่าง xmin และ xmax def bisection4(f,xmin, xmax): "พยายามหาค่า x ที่ทำให้ f(x)==0, xmin และ xmax คือช่วงที่เดาว่าคำตอบอยู่ในนั้น" tolerance = 1e-6 #ตั้งค่า tolerance ไว้ให้หยุดทำงาน #จะหยุดทำงานเมื่อ xmin ห่างจาก xmax น้อยกว่า tolerance if xmin > xmax: # จัดการให้ xmin น้อยกว่า xmax เสมอ xmin, xmax = xmax, xmin if f(xmin)*f(xmax) > 0: print("Bad xmin, xmax") #ถ้าฟังก์ชั่น f ไม่เปลี่ยนเครื่องหมายจากบวกไปลบหรือลบไปบวก #ระหว่าง xmin และ xmax ก็ควรไปเดาใหม่ว่า xmin, xmax คืออะไร return None while xmax-xmin > tolerance: # ทำตรงนี้วนๆไปตราบใดท่ี xmin และ xmax ยังห่างกันกว่า tolerance #print(f"xmin, xmax = {xmin}, {xmax}\t f(xmin)*f(xmax) = {f(xmin)*f(xmax):.5f}") if f(xmin) == 0: return xmin if f(xmax) == 0: return xmax xmid = (xmin + xmax)/2 #หาจุดกลางระหว่าง xmin และ xmax if f(xmin)*f(xmid) < 0: #ถ้าฟังก์ชั่นเปลี่ยนเครื่องหมายระหว่าง xmin กับจุดกลาง ก็ไปหาต่อในช่วงนี้ xmin, xmax = xmin, xmid else: #ไม่งั้นก็หาต่อในช่วงจุดกลางถึง xmax xmin, xmax = xmid, xmax slope = (f(xmax)-f(xmin))/(xmax-xmin) x = xmin + (0-f(xmin))/slope # คำตอบที่ได้จะเทียบบัญญัติไตรยางค์ระหว่างจุด (xmin, f(xmin)) และ จุด (xmax, f(xmax)) return x # In[14]: # ทดลองเดา xmin, xmax = -3, -2 x = bisection4(quintic,-3,-2) print(f"Value at x = {x:.6f} is {quintic(x):.6f}") # In[15]: # ทดลองเดา xmin, xmax = -3, 0 x = bisection4(quintic,-3,0) print(f"Value at x = {x:.6f} is {quintic(x):.6f}") # ลองหาคำตอบ $x^2 - 4 = 0$ แถวๆระหว่าง 1.7 และ 2.2 (คำตอบที่ถูกคือ 2) # ![Screen%20Shot%202562-06-29%20at%2016.32.42.png](attachment:Screen%20Shot%202562-06-29%20at%2016.32.42.png) # In[16]: # ลองใช้วิธี bisection หาคำตอบของควอดราติก x**2 - 4 = 0 แถวๆ (1.7,2.2): def quad(x): return x**2 - 4 x = bisection4(quad, 1.7, 2.2) print(f"Value at x = {x:.6f} is {quad(x):.6f}") # ลองหาคำตอบ $\frac{1}{1+x^2} - 1/2 = 0$ แถวๆระหว่าง 0.7 และ 1.2 (คำตอบที่ถูกต้องคือ x = 1) # # ![Screen%20Shot%202562-06-29%20at%2016.12.32.png](attachment:Screen%20Shot%202562-06-29%20at%2016.12.32.png) # In[17]: # หาคำตอบ 1/(1+x**2) - 1/2 = 0 แถวๆระหว่าง 0.7 และ 1.2 (คำตอบที่ถูกต้องคือ x = 1) def f1(x): return 1/(1+x**2) - 1/2 x = bisection4(f1, 0.7, 1.2) print(f"Value at x = {x:.6f} is {f1(x):.6f}") # ลองหาคำตอบ $cos(\frac{x}{2}) = 0$ แถวๆ 2.8 และ 3.2 (คำตอบที่ถูกต้องคือค่า $\pi$) # ![Screen%20Shot%202562-06-29%20at%2016.18.08.png](attachment:Screen%20Shot%202562-06-29%20at%2016.18.08.png) # In[18]: def f2(x): return math.cos(x/2) x = bisection4(f2, 2.8, 3.2) print(f"Value at x = {x:.6f} is {f2(x):.6f}") # In[ ]: