มีเด็กๆที่สนใจคณิตศาสตร์มาถามผมว่าสมการทั่วๆไปแก้ออกมาเป๊ะๆไม่ได้แล้วเราทำอย่างไร ผมก็บอกว่าสมการส่วนใหญ่เราต้องหาคำตอบด้วยการประมาณเอาครับ ซึ่งวิธีอันหนึ่งที่เราสามารถใช้ได้ง่ายๆก็คือวิธีของนิวตัน ข้างล่างนี้เป็นกระทู้ที่ผมเขียนไว้ที่ Mahidol Physics Educational Center ครับ:
สำหรับปัญหาที่เราแก้สมการโดยตรงไม่ได้ เราต้องแก้ด้วยวิธีประมาณด้วยตัวเลขครับ วิธีที่ใช้กันบ่อยๆวิธีหนึ่งก็คือวิธีการของนิวตัน (Newton’s method: http://en.wikipedia.org/wiki/Newton’s_method) ครับ
วิธีการของนิวตันบอกว่า ถ้าจะแก้สมการ f(x) = 0 ให้เราเดาค่า x มาสักค่า (เรียกมันว่า x0) ก็แล้วกัน แล้วค่า x อันต่อไป (เรียกมันว่า x1) ที่น่าจะทำให้ f(x) ใกล้ศูนย์มากขึ้น ควรจะคำนวณอย่างนี้ครับ:
x1 = x0 – f(x0)/f'(x0) โดยที่ f'(x) คือ derivative ของ f(x) ครับ
ถ้าค่า x1 ทำให้ f(x) ไม่ใกล้ 0 พอ เราก็หา x2, x3, x4, … ไปเรื่อยๆจนเราพอใจว่าค่า f(xn) ใกล้ 0 พอแล้ว โดยที่ xn หาได้จาก xn-1 ดังนี้ครับ:
xn = xn-1 – f(xn-1)/f'(xn-1)
ถ้าจะทำการคำนวณด้วยวิธีของนิวตันใน Mathematica สามารถทำอย่างนี้ครับ:
newtonSolve[f_, guess_, steps_] := NestList[ #1 – f[#1]/f'[#1] &, guess, steps]
f คือฟังค์ชั่นที่เราจะหา f(x) = 0
guess คือค่าที่เราเดาตอนแรกว่า f(guess) น่าจะไม่ห่างจาก 0 นัก
steps คือจำนวนครั้งที่เราจะทำการทำวิธีของนิวตันซำ้ๆกัน
ยกตัวอย่างเช่น เราจะหาค่ารูทที่สองของสอง เราก็เขียนสมการ f(x) = x^2 -2 = 0 ก่อน เพราะค่า x เท่ากับรูทที่สองของสองจะแก้สมการนั้นพอดี:
f[x_] := x^2 – 2
แล้วเราก็เรียก newtonSolve ด้วยฟังค์ชั่น f โดยเดาค่า guess = 1 และให้ทำซ้ำสักห้าครั้ง:
newtonSolve[f, 1, 5]
แล้วเราก็ได้ผลดังนี้: {1, 3/2, 17/12, 577/408, 665857/470832, 886731088897/627013566048}
Mathematica ทำการคำนวณให้เป็นค่าเศษส่วนไม่มีทศนิยม เพราะเราเดาด้วยค่า 1 ซึ่งเป็นจำนวนที่ไม่มีการประมาณเข้ามาเกี่ยวข้อง ถ้าเราต้องการคำตอบเป็นเลขทศนิยม เราก็ควรเดาด้วยค่า 1.0 ดังนี้:
newtonSolve[f, 1.0, 5]
แล้วเราก็จะได้ผลดังนี้: {1., 1.5, 1.41667, 1.41422, 1.41421, 1.41421} ซึ่งเราจะเห็นว่า 1.41421 นั้นเป็นค่าประมาณของรูทที่สองของสองได้ดีทีเดียว
นักศึกษาลองไปทดลองดูครับ
มีใครอยากลองอธิบายว่า newtonSolve[f_, guess_, steps_] := NestList[ #1 – f[#1]/f'[#1] &, guess, steps] ทำงานอย่างไรใน Mathematica ไหมครับ เป็นการฝึกความเข้าใจเรื่อง Pure function และ Functional programming ครับ
