บทที่ 6, 7 พหุนาม และ อินเตอร์โพเลชั่น

Chapter 6

พหุนาม

6.1 ส่วนนำเรื่อง

            เอกสารนี้ออกแบบเพื่อให้ผู้อ่านคุ้นเคยกับการหาค่าพหุนามและการอินเตอร์โพเลชัน (interpolation) และการใช้  least squares fitting ของข้อมูลกับพหุนาม.

6.2 การหาค่าพหุนาม

            ใช้ไซแลบในการคำนวณฟังก์ชันต่อไปนี้ได้อย่างชัดเจนปราศจากจัดเรียกให้หรือการทำ

แฟคเตอร์

f(x) = x⁸ − 8 x⁷ + 28 x⁶ − 56 x⁵ + 70 x⁴ − 56 x³ + 28 x² − 8^x + 1

g(x) = (((((((x − 8) x + 28) x − 56) x + 70) x − 56) x + 28) x − 8) x + 1

h(x) = (x − 1)⁸

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

แนะ: อาจต้องการสร้างไฟล์หนึ่งที่ฟังก์ชันทั้งสาม f,g และ h กำหนดขึ้นในไซแลบ  ในไซแลบคำสั่ง x^4 จะพยายามที่จะยกระดับเมทริกซ์ x  เข้าสู่ยกกำลัง 4  (การใช้การคูณเมทริกซ์)  ในอีกทางหนึ่งตัวกระทำด็อท x.^4 จะใช้กับตัวกระทำ ^ ไปยังการนำเข้าเมทริกซ์ x  ทำให้แนะนำได้ว่า ผู้ใช้มักจะใช้เครื่องหมายด็อท เมื่อกำหนดฟังก์ชันดังกล่าว.

6.3 พหุนาม

            การทำให้ขั้นตอนวิธีของไซแลบจัดการกับพหุนาม
พหุนาม a1 +a2x + · · · + amx^m−1 สามารถที่จะแทนในไซแลบตามสัญลักษณ์ดังนี้
ในตัวอย่างต่อไป -4 - 3x + x² แทนโดยเว็คเตอร์สัมประสิทธิ์  [ -4 - 3  1]  โดยพหุนาม p ที่กำหนดขึ้นโดยคำสั่งดังนี้
v = [-4,-3,1]
p = poly(v,’x’,’coeff’)

            การสร้างพหุนามสามารถทำได้ในแนวทางที่เป็นธรรมชาติ ดังเช่น

x = poly(0,’x’)               //  แก่นเมล็ดของพหุนามที่ใช้ตัวแปร ’x’

p = x^2 - 3*x - 4             // แทนการนำเสนอ x^2 - 3 x -4

มีฟังก์ชันมากมายที่มีประโยชน์ ที่สามารถประยุกต์ใช้กับพหุนาม สามารถหารากเป็นจำนวนเชิงตัวเลขของพหุนาม  ให้หาค่าพหุนามหนึ่ง หรือ ค่าค่าอนุพันธ์

z = roots( p )

z(1)^2 -3*z(1) -4            //  สองแนวทางในการหค่าโพลิโนเมียล

horner(p,z(1))                // ค่าแรกที่ 1

derivat(p)                      // คำนวณอนุพันธ์ของโพลีโนเมียล

            วิธีการโพลีโนเมียลของไซแลบได้ให้ฟังก์ชันที่เป็นสัดส่วน ตัวอย่างเช่น

x = poly(0,’x’)               // แก่นของพหุนามใช้ตัวแปร ’x’

p = x^2 - 3*x - 4             // แทน x^2 - 3 x –4

r = x/p                           //แทนสัดส่วน

                                    // ฟังก์ชัน x/(x^2 - 3 x -4)

horner(r,1.0)                   // หาค่า r ที่ x=1

6.4 Linear Least Squares (Heath Computer Problem 3.6)

            ส่วนของคำถามเป็นส่วนของแบบฝึกหัดที่สอง (of assignment 2) เรื่องนี้ให้โอกาสในการทดลองหาว่าปัญหาด้วยเอกสารช่วยเหลือนี้จะช่วยให้อยู่เหนือการ “bumps”.

            เพื่อที่จะสาธิตความแตกต่างของจำนวนระหว่างวิธีสมการปกติกับQR factorisation สำหรับ  least squares เชิงเส้น  นั้นต้องการปํยหาที่เงื่อนไขคลุมเคลือ และมักจะมีส่วนเหลือเล็กน้อย  นั้นสามารถสร้างปัญหาเช่นนั้นให้เป็นดังต่อไป คือ  จะได้ฟิตโพลีโนเมียลที่ degree n -1

P_{n -1}(t) = x1 + x2t + x3t² + · · · + xnt^n-1

            สำหรับ  m จุดข้อมูล (ti, yi), ที่ m > n.  เลือกให้ ti = (i - 1)/(m - 1), i = 1, ...,m,

ดังนั้นจุดข้อมูลมีช่องว่างเท่ากันในระหว่างช่วงข้อมูลl [0, 1].

            สร้างให้ข้อมูลบางอย่าง เป็น y   ค่าของ yi  สร้างขึ้นจากการเลือกค่าแรกของ  กล่าวคือ xj = 1, j = 1, ..., n, และค่าค่าผลลัพธ์ของโพลีโนเมียลให้ได้ค่า yi = p_n-1(ti), i = 1, ...,m.  โดยต้องการจะดูว่าถ้าสามารถคืนค่าxj ที่ใช้ในการสร้าง yi, แต่เพื่อทำให้น่าสนใจยิ่งขึ้น  แรกสุดกระตุ้นให้เกิดค่า yi แบบสุ่ม จำลองแบบข้อมูลที่ผิดพลาดซึ่งเป็นปัญหาตามปกติของลีสท์สแคว์  โดยเฉพาะอย่างยิ่ง ที่นำเอา  yi = p_n-1(ti) +

(2u_i − 1) ∗ , i = 1, ...,m,  เมื่อ u_i คือจำนวนตัวเลขสุ่มที่แจกแจงอย่างสม่ำเสมอในช่วง

[0, 1] และ  เป็นจำนวณค่าน้อยที่ที่หาค่าการรบกวนสูงสุด ด้วยค่าความละเอียดสองเท่าของ IEEE  พารามีเตอร์ที่มีเหตุผลลสำหรับปัญหาคือ  m = 21,n = 12, and  = 10−10.

            ต่อไปจะได้เปรียบเทียบ 2 แบบในการคำนวณผลเฉลย least square  สำหรับปัญหาการฟิตข้อมูลโพลีโนเมียล

1. สร้างชุดข้อมูล (ti, yi), เป็นเหมือนหัวรายการ ให้ใช้ฟังก์ชัน rand เพื่อสร้างข้อผิดพลาดสุ่ม

2. สร้างเมทริกซ์ Vandermonde  A สำหรับค่าทั่วไป m และ n.

3. สร้างเว็คเตอร์ข้อมูลทางด้านขวามือ b  ประกอบด้วยค่าข้อมูล y

4. ให้แก้ปัญหาระบบ least squares โดยใช้สมการตามปกติสำหรับปัญหานี้

      เช่นแก้สมการ (ATA)x = ATb for x.

5. การแก้ปัญหาระบบ least squares โดยใช้วิธี orthogonal factorisation

(Form the QR factorisation using the Scilab command qr.)

6. คำสั่ง backslash ของไซแลบสามารถใช้สำหรับระบบที่หาได้เกินจำเป็น Ax ≈ b. ให้ใช้วิธีการนี้ และเปรียบเทียบผลกับวิธีสมการปกติ และวิธี QR  วิธีการใดที่คิดว่าคำสั่ง Scilab backslash ใช้สำหรับระบบที่เกินจากที่หาได้

7. จงเปรียบเทียบผลเฉลย 3 อย่างในเว็คเตอร์ วิธีการใดให้ผลเฉลยที่ไวต่อการรบกวนที่นำเข้ามาในข้อมูล? วิธีการใดเข้าใกล้การคืนกลับค่า x ที่ใช้ในการสร้างข้อมูล?

8. โดยการวัดค่าขนาดของข้อผิดพลาดในผลเฉลย  ที่ให้ขอบเขตล่างสำหรับเงื่อนไขสำหรับแต่ละวิธีการ.

9. ให้ค้นคว้าสืบเสาะดูว่าความแตกต่างระหว่างผลเฉลยผลต่อความสามารถในการลงจุด (ti,yi) ท่เหมาะสมใกล้เคียงกับพหนุนามหรือไม่  เพราะอะไร?

Chapter 7

Interpolation

7.1 นำเรื่อง

            เอกสารเสริมนี้ออกแบบเพื่อให้ผู้อ่านได้คุ้นเคยกับการใช้  interpolation.

7.2 Interpolating Runge’s Function

            พิจารณาฟังก์ชัน

f(t) =  1/ 1 + 25t² .

            ฟังก์ชันนี้รู้จักกันในนามฟังก์ชัน  Runge’.

            กำหนดฟังก์ชันไซแลบเพื่อคำนวณฟังก์ชนของ Runge’

function y=runge(t)

//

// Runge’s Function

// Standard example of polynomial interpolation

// creating oscillations

//

y=(1.0)./(1+25*t.^2)

endfunction

ข้อสังเกตที่นำเอาวงเล็บรอบ (1.0)  เพื่อให้แน่ใจว่านำคำสั่ง ./ มาใช้  อันเป็นนิพจน์ที่อยู่ในรูปแบบ1./(..)  จริงแล้วคือการใช้ตัวกระทำเมทริกซ์ /

            พิจารณาพหุนาม

            p_n(t) =

            ซึ่งการประมาณค่าในช่วง (ค่าที่ตรงกับ)ฟังก์ชัน Runge ที่ n ที่จุดที่ว่างเท่ากัน ระหว่างค่า-1 และ 1 (นั่นคือที่จุด -1:2/(n-1):1|.   โดยมีตัวไม่ทราบค่า  n จำนวน xi  สำหรับค่าสัมประสิทธิ์ของพหุนามและจำนวน n จุดข้อมูล  ที่ให้สมการเชิงเส้นสำหรับสัมประสิทธิ์พหุนามที่ได้ให้ค่าจุดข้อมูล   ให้สร้างฟังก์ชันหนึ่งเรียกว่า polyfit, ด้วยกากรเรียกใช้ sequence pp=polyfit(tdata,ydata)  ที่สอดรับกับพหุนามดีกรี n - 1 ถึง n จุดข้อมูล t, y   และคืนค่าพหุนาม pp.   สามารถที่จะเพิ่มอาร์กิวเมนต์ พิเศษได้ง่ายในการคำนวณ least

square fit พหุนามดีกรี m < n ไปยัง n จุดข้อมูล.  ผลลัพธ์ของ polyfit นั้นแล้วสามารถที่จะหาค่าโดยใช้คำสั่งดังนี้

tt = -1:0.01:1;

yy = horner(pp,tt)

แนะ:  สามารถทำได้โดยใช้ตัวกระทำ แบคสแลทมาตรฐานเพื่อแก้ปัญหาสัมประสิทธิ์  ทันทีที่มีเมทริกซ์ Vandermonde  คำสั่ง feval กับการกำหนดฟังก์ชันที่เหมาะสมกดังนี้

function z=vandermonde(t,q)

z=t.^q;

endfunction

            อาจจะมีประโยชนในการสร้างเมทริกซ์ Vandermonde ที่เกี่ยวพันธ์กัน  โดยอาจจะชอบที่จะใช้deff เพื่อกำหนดฟังก์ชันขี้นเอง

deff(’z=vandermonde(t,q)’,’z=t.^q’)

            ข้างบนนี้มีประโยชน์สำหรับการใช้ฟังก์ชันสั้นๆ   การใช้ฟังก์ชัน polyfit ในการคำนวณสัมประสิทธิ์พหุนามถึงดีกรี 5 และ ดีกรี 10 ซึ่งปรมาณค่าช่วงฟังก์ชัน Runges ที่ 6 และ 11 ของจุดห่างเท่ากันในช่วง [−1, 1]. การประมาณการของพหุนามประสบผลสำเร็จอย่างดี เกี่ยวกับศูนย์และจุดปลาย
7.3  พหุนามเทเลอร์ (Taylor Polynomials)

            ข้อสังเกตที่ว่านิพจน์ 1/(1 + t) = 1 − t + t2 − t3 ง่ายต่อการให้พหุนามเทเลอร์มีศูนย์กลางโดยรอบจุดกำเนินของฟังก์ชัน Runge จัดอยู่ในรูปแบบคือ

            1 − 25t2 + (25t2)2 − (25t2)3...

            ให้ใช้รูปแบบดังกล่าวเพื่อพล็อตพหุนามเทเลอร์องศาที่ 4 โดยมีศูนย์กลางที่ศูนย์สำหรับฟังก์ชัน Runge คลุมอยู่ระหว่างช่วง [-1,1].  ควรจะได้พล็อตฟังก์ชัน Runge และพหุนามเทเลอร์คลุมช่วงที่เล็กลงกล่าวคือ [-0.2,0.2]. ให้พยามอธิบายถึงค่าความละเอียดถูกต้อง และอะไรที่คาดหวังถึงความละเอียดถูกต้องที่องศาที่ 8 ของพหุนามว่าจะเหมือนอะไร

7.4 การคาดคะเนช่วงแบบ Chebyshev

            ให้มาดูการคำนวณและการประมาณพหุนามกราฟ กับฟังก์ชัน Runge ไปยังองศาที่ 4 และ 8 โดยใช้การคาดคะเนช่วง ที่จุดChebyshev  สำหรับกรณี องศาที่ 4 นิพจน์ไซแลบที่ใช้คำนวณ 5 จุดที่จำเป็นคือ

data = cos((1:2:(2*5-1))*%pi/(2*5))

ydata = runge(tdata)

chebp = polyfit(tdata,ydata)

tt = -1:.01:1;

rr = runge(tt);

yy = horner(chebp,tt);

xbasc();

plot(tt,yy,style=1);

plot2d(tt,rr,style=2);

plot2d(tdata,ydata,style=-1)

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

            คำสั่งไซแลบต่อไปใช้คำสั่ง splin เพื่อคำนวณการประมาณค่าช่วงผ่านจุดที่ระบุโดยเว็คเตอร์ x และ y.  โดยใช้คำสั่ง interp เพื่อหาค่า spline ที่จุด xx.

tdata = -1:.5:1;

ydata = runge(xs);

ddata = splin(tdata,ydata);

tt = -1:.01:1;

rr = runge(tt);

ss = interp(tt,tdata,ydata,ddata);

Thus the spline approximation can then be graphed using

xbasc();

plot2d(tt,ss,style=1)

plot2d(tt,rr,style=2);

plot2d(tdata,ydata,style=-1);

            ให้เปลี่ยนรหัสคำสั่งในการประมาณการพล็อต cubic splineของฟังก์ชัน  Runge  โดยใช้ชุดของจุด9 จุดที่ช่วงว่างห่างเท่าๆ กันในการประมาณค่าช่วงของฟังก์ชัน Runge

7.6 ประมาณค่าช่วงแบบ Monotonic Cubic

From the previous example we see that cubic spline interpolation can still in-

troduce unwanted oscillations. There are many techniques, for instance a fairly

common method is that of Fritsch and Carlson (SIAM J. Numer. Anal. 17,

pp 238-246 1980). Scilab implements a monotonic interpolation method as an

option to the splin function.

Here is some data

tdata = 0:1:10;

ydata = [0 0 0 0 0.1 1.8 1.9 2.0 2.0 2.0 2.1];

Lab Book: Fit an ordinary spline function to this data set. Produce a plot

I obtained the following:

0 1 2 3 4 5 6 7 8 9 10

-0.2

0.2

0.6

1.0

1.4

1.8

2.2

+ + + +

+

+

+

+ + +

+

The result is quite oscillatory, maybe not really what we want.

Now re do the spline interpolation of this data, but use the extra argument

’monotone’. I.e. use

ddata = splin(tdata,ydata, ’monotone’);

Now you can use tdata and ydata values together with the new derivative

values to produce an interpolant using the interp function.

Lab Book: Plot the data points, together with plots of the ordinary spline and

monocubic interpolant. Comment on the respective interpolants.

Hopefully you can see that in some cases it is useful to try to maintain the

shape implicit in your data sets. This is particularly true for coarse data sets.