Khớp đường cong và thuật toán nội suy Spline
30/09/2008

Thuật ngữ đường cong (curvefitting) hay điều chỉnh dữ liệu được dùng để mô tả bài toán tổng quát của việc tìm các hàm khớp với một tập các giá trị đang được quan sát ứng với một tập điểm.


Cho các điểm:

x1, x2,…, xN

và các giá trị tương ứng

y1, y2,…, yN

Mục đích là tìm các hàm sao cho

f(x1) = y1, f(x2)=y2, …, f(xN)=yN

và sao cho f(x) được giả sử “hợp lý” ở những điểm dữ liệu khác. Có thể rằng một x và y được liên hệ bởi hàm f(x) ẩn danh nào đó và mục tiêu của ta là tìm ra hàm đó, nhưng định nghĩa của từ “hợp lý” phụ thuộc vào từng ứng dụng. Ta thấy rằng thường thì xác định hàm “không hợp lý” thì dễ hơn.

Khớp đường cong có ứng dụng hiển nhiên trong sự phân tích các dữ liệu thuộc thí nghiệm và còn nhiều ứng dụng khác nữa…Ví dụ, nó có thể dùng trong đồ hoạ máy tính để sản sinh ra đường cong “coi được” mà không cần pảhi lưu một số lượng lớn các điểm vẽ. Một ứng dụng có liên hệ là dùng chỉnh đường cong để cho ra một thuật giải nhanh trong tính toán giá trị của hàm chưa biết ở một điểm bất kỳ: Giữ một bảng nhỏ chứa các giá trị chính xác, sự hiệu chỉnh đường cong sẽ suy ra các điểm khác.

Bài viết này sẽ đề cập tới một trong các phương pháp cơ bản để tiếp cận bài toán này đó là phương pháp nội suy: tìm một hàm liên tục khớp với các giá trị đã cho.

Các đa thức bậc thấp là những đường cong đơn giản được sử dụng rộng rãi trong nối đường cong. Thay vì dùng các đa thức khác nhau để nối các điểm kề nhau, nối các đoạn sao cho thật mịn. Một trường hợp đặc biệt liên hệ sự tính toán tương đối trực tiếp, phương pháp này còn gọi là nội suy spline.

Spline là một thiết bị cơ học được người vẽ sơ đồ thiết kế dùng để vẽ các đường cong đẹp, có thẩm mỹ: chỉ cần xác định tập hợp các điểm (nút) rồi bẻ cong một giải plastic hay miếng gỗ linh hoạt (spline) quanh chúng và lấy vết chúng để tạo thành một đường cong. Nội suy spline thì tương đương về mặt toán học với tiến trình này và cho ra cùng một kết quả. Hình dưới minh hoạ một spline qua 10 điểm. 

 

Có thể thấy rằng hình dạng của một đường cong tạo bởi spline giữa hai nút kề nhau là một đa thức bậc ba. Trở lại bài toán nối dữ liệu, điều này có nghĩa là ta nên xem đường cong là N-1 đa thức khác nhau có bậc ba.

si(x) = aix3 + bix2+ cix + di i = 1, 2…N-1

Với si(x) là đa thức bậc ba xác định giữa khoảng xi và xi+1. Spline có thể được biểu diễn trong một mảng bốn chiều (hay trong 4x(N-1) mảng hai chiều). Việc tạo một spline gồm việc tính các hệ số a, b,c, d từ các điểm x và các giá trị y đã cho. Việc bẻ cong một spline về mặt vật lý tương ứng với việc giải hệ phương trình với nghiệm là các hệ số. Ví dụ, Hiển nhiên ta phải có: si(xi) = yi và si+1(xi+1) = yi+1 với i = 1, 2, …, N-1 vì spline phải chạm vác nút. Không những nó phải chạm mặt toán học được này nghĩa là các đạo hàm bậc nhất của các đa thức spline phải bằng nhau ở mỗi nút (s’i-1(xi)= s’i(xi) với i=2, 3, …N-1). Thật sự thì các đạo hàm bậc hai của các đa thức cũng phải bằng nhau ở mỗi nút. Các điều kiện này cho ra 4N-6 phương trình với 4(N-1) hệ số là ẩn. Cần xác định thêm hai điều kiện nữa để mô tả tình trạng ở hai điểm cuối của spline. Có nhiều cách: dùng cái gọi là spline ‘tự nhiên’ được rút từ s”i(xi) = 0 và s”N-1(xN) = 0. Các điều kiện này cho ra một hệ đầy đủ 4N-4 phương trình với 4N- 4 ẩn số, hệ phương trình này giải được bằng phép khử Gauss với ẩn là các hệ số.

Tuy nhiên, cũng cùng một spline nhưng có thể tính toán khá hiệu quả hơn vì thật sự chỉ N-2 ẩn: hầu hết các điều kiện của spline là thừa. Ví dụ, giả sử pi là đạo hàm bậc hai của spline tại điểm xi, vì s”i-1(xi) = s”i(xi) = pi với i=2, 3, …, N-1 và pi= pN = 0 nếu biết trước các giá trị pi và pN, thì tất cả các hệ số a, b, c, d có thể được tính trên các đoạn spline, vì ta có bốn phương trình với bốn ẩn số trên các đoạn: với i=1, 2, …,N-1 ta phải có:

si(xi) = yi ; si(xi+1) = yi+1 ; s”i(xi) = pi ; s”i(xi+1) = pi+1

Các giá trị x, y đã cho trước, để xác định đầy đủ spline chỉ cần tính các giá trị p2, …, pN. Để tính được, dùng điều kiện đạo hàm bậc nhất phải bằng nhau: có N-2 điều kiện ứng với N-2 phương trình cần để giải N-2 ẩn, theo pi .Để diễn tả các hệ số các hệ số a, b, c, d theo p, các giá trị đạo hàm bậc hai, ta thay các biểu thức vào bốn phương trình trên trên mỗi đoạn spline, điều này dẫn đến một số biểu thức phức tạp không cần thiết. Thay vì diễn tả các phương trình trên từng đoạn spline ở dạng chuẩn liên hệ đến ít ẩn số hơn. Nếu ta thay các biến là t = (x-xi)/(xi+1-xi) thì spline được biểu diễn như sau:

si(t) = tyi+1 + (1-t)yi + (xi+1 -xi)2((t3- t)pi+1 – ((1-t)3 – (1-t)pi)/6

bây giờ mỗi spline được xác định trên đoạn [0,1]. Phương trình này đơn giản hơn là hình thức bên ngoài của nó, vì ta chỉ quan tâm đến hai điểm nút 0 và 1 và hoặc t hoặc (1-t) bằng 0 ở những điểm mút này. Cách biểu diễn này dễ kiểm spline được nội suy và liên tục vì si-1(1) = si(0) = yi với i = 2,…, N-1, nó chỉ hơi khó chứng minh bậc hai liên tục vì s”i(1) = s”i+1(0)= pi+1. Có các đa thức bậc ba thoả những điều kiện này tại các điểm mút, vì thế chúng tương đương với các đoạn spline mô tả ở trên. Nếu ta thay thành t và tìm hệ số của x3,… thì ta sẽ có cùng biểu thức ẩn là a, b, c, d theo x, y và p giống như ta dùng phương pháp được mô tả ở đoạn trước. Nhưng không có lý do gì làm như vậy, vì ta đã kiểm rằng các đoạn spline này thoả các điều kiện cuối, và có thể lượng giá từng spline ở bất kỳ điểm nào trong đoạn của nó bằng cách tính t và dùng công thức trên (khi đã biết p).

Để tính p, ta gán các đạo hàm cấp một của các đoạn spline bằng các điểm cuối. Đạo hàm cấp 1 (theo x) của các phương trình trên là:

s’i(t) = zi + (xi+1 - xi)((3t2 - 1)pi+1– (3(1-t)2 – 1)pi)/6

Trong đó zi= (yi+1 - yi)/(xi+1 - xi). Kế đến gán s’i-1(1) = s’i(0) với i = 2…N-1, ta sẽ có hệ N-2 phương trình:

(xi-xi-1)pi-1 + 2(xi+1 – xi-1)pi + (xi+1 - xi)pi+1 = 6(zi – zi-1)

Hệ phương trình này thuộc dạng tam giác đơn giản dùng phép khử Gauss để giải. Ví dụ, nếu đặt ui = xi+1– xi, di = 2(xi+1 – xi-1), và wi = 6(zi – zi+1), ta có hệ phương trình với N=7:

 

Thật sự, hệ tam giác đối xứng này có nửa dưới bằng với nửa trên của đường chéo chính. Nó được xoay trên phần tử có sẵn lớn nhất nếu khong cần tìm lời giải đúng cho hệ phương trình này.

Phương pháp được mô tả ở đoạn t rên để tính spline bậc ba được dễ dàng viết dưới dạng Pascal như thủ tục makesplinedưới đây:

++++

procedure makespline;

var i: integer;

begin

readln(N);

for i:= 1 to N do readln(x[i], y[i]);

for i:= 2 to N-1 do d[i]: = 2 + (x[i+1] - x[i-1]);

for i:= 1 to N-1 do u[i]: = x[i+1] - x[i];

for i:= 2 to N-1 do

w[i]: = 6.0 * ((y[i+1) - y[i])/u[i] – (y[i] - y[i-1])/u[i-1]);

p[1]:= 0.0; p[N]:= 0.0;

for i:= 2 to N-2 do

begin

w[i+1]:= w[i+1] - w[i]* u[i]/d[i];

d[i+1]:= d[i+1] - u[i]* u[i]/d[i];

end;

for i:= N-1 downto N-1 do

p[i]:= (w[i] - u[i]*p[i+1])/d[i];

end;

+++

Các mảng d và u là đại diện matrận tam giác.

*** Mộ spline bậc ba trên N điểm có thể thực hiện theo thời gian tuyến tính.

Điều này hiển nhiên có được từ chương trình.

Ví dụ xây dựng một spline bậc ba từ năm điểm: (1.0, 2.0), (2.0, 1.5), (4.0, 1.25), (5.0,1.2), (8.0, 1.125), (10.0, 1.1).

(Các điểm trên lấy từ hàm 1+1/x). Các tham số spline tìm được bằng cách giải hệ phương trình:

 

với kết quả: p2 = 0.39541, p3 = -0.06123, p4 = 0.02658, p5 = -0.00047

Để lượng giá spline cho bất kỳ giá trị x nào trong khoảng [x1, xN], đơn giản ta tìm khoảng [xi, xi+i] có chứa x, rồi tính t và dùng công thức trên cho si(x) (rồi dùng lại hàm này để tính pi và pi+1).

 

function eval(v: real):real;

var t: real; i: integer;

funtion f(x: real):real;

begin f:= x*x*x-x; end;

begin i:= 0;

repeat i:=i+1 until v[i+1];

t:= (u-x[i])/u[i];

eval:=t*y[i+1]+(1t)*y[i]+u[i]*(f(t)*p[i+1] + f(t-1]+f(t-1)*p[i] )/6.0;

end;

Chương trình này không kiểm tra khi v không nằm giữa x[1] và x[N]. Nếu số đoạn spline nhiều (nghĩa là N lớn) có một số phương pháp tìm có hiệu quả hơn.

Có nhiều thay đổi tư tưởng điều chỉnh đường cong bằng cách ráp các đoạn đa thức sao cho thật mịn. Việc tính toán các spline thuộc lĩnh vực nghiên cứu đã phát triển rất tốt. Các kiểu spline khác liên hệ đến các chuẩn làm mịn cũng như các thay đổi khác như nới lỏng điều kiện các spline phải chạm đúng vào mỗi điểm dữ liệu. Về mặt tính toán, chúng liên hệ cùng một số bước xác định các hệ số cho mỗi đoạn spline bằng cách giải hệ phương trình tuyến tính rút từ các giới hạn được đưa ra trên sự kiện chúng được nối với nhau như thế nào.
URL của bài viết này::http://www.vnschool.net/modules.php?name=News&file=article&sid=2502