Hướng dẫn Vẽ đường cong và Nội suy bằng Matlab phần 2

1. Khớp đường cong:

Mời các bạn xem phần 1: Hướng dẫn Vẽ đường cong và Nội suy bằng Matlab phần 1

Ở bài trước chúng ta đã hiểu thế nào là khớp đường cong thông qua đoạn code sau:

>> x = [0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1];
>> y = [-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
>> n = 2;
>> p = polyfit(x,y,n);
>> p

p =

   -9.8108   20.1293   -0.0317

>> xi = linspace(0,1,100);
>> z=polyval(p,xi);
>> plot(x,y,'o',x,y,xi,z,':')
>> xlabel('x'),ylabel('y = f(x)')
>> title('Xap xi duong cong bac 2')

2. Nội suy bậc 1.

Như đã nói trong mục 1, nội suy là cách đánh giá hàm tại các điểm nằm giữa các điểm số liệu. Nội suy là một công cụ quý giá khi ta không thể nhanh chóng đánh giá hàm tại các điểm trung gian mong muốn. Chẳng hạn điều này rất đúng cho trường hợp các điểm số liệu là kết quả đo đạc thí nghiệm hoặc thủ tục tính toán dài dòng.

Theo ngầm định MatLab nối các điểm dữ liệu bằng đường thẳng để tạo nên đồ thị. Nội suy tuyến tính giả thiết là các giá trị trung gian nằm trên đường thẳng nối các điểm được nhập vào. Chắc chắn là khi số các điểm dữ liệu tăng lên và khoảng cách giữa chúng giảm đi thì nội suy tuyến tính trở nên chính xác hơn. Ví dụ :

>> x1=linspace(0,2*pi,60);
>> x2=linspace(0,2*pi,6);
>> plot(x1,sin(x1), x2,sin(x2),'-')
>> xlabel('x'), ylabel('y=sin(x)')
>> title(' Noi suy tuyến tính")

Trong hai đồ thị của hàm sin nhận thấy rõ ràng là tại các vị trí nằm giữa các điểm số liệu, thì đường dùng 60 điểm dữ liệu sẽ trơn nhẵn và chính xác hơn nhiều so với đường dùng 6 điểm dữ liệu.
Cũng như với khớp đường cong, ở đây ta cũng phải lựa chọn. Có nhiều cách tiếp cận để nội suy dựa trên các giả thiết khác nhau. Hơn nữa lại có khả năng nội suy không chỉ một chiều. Tức là nếu bạn có dữ liệu phản ánh hàm 2 biến z = f(x,y), bạn có thể nội suy giữa các giá trị của x và y để xác định giá trị trung gian của z. MatLab cung cấp cho bạn cách tuỳ chọn nội suy bằng hàm một chiều interp1 và hàm hai chiều interp2.
Để minh hoạ nội suy một chiều, ta xét trường hợp sau.
Thực hiện một phần trong đề tài khoa học ghi nhiệt độ ở Trường Xuân từng giờ một trong thời gian 12 tiếng để sử dụng thông tin này trong báo cáo về tình hình khí hậu địa phương. Phân tích dữ liệu :

>> hours = 1:12;
>> temps=[5 8 9 15 25 29 31 30 22 25 27 24] %là các lần đọc nhiệt độ (đo bằng °C).

Dựa vào các số liệu này, có thể vẽ đồ thị như sau :

>> plot(hours, temps, hours, temps, '+')
>> xlabel('hours'), ylabel(' temps')
>> title(' Nhiet do o Xuan Truong')

Nhiệt độ ở Trường xuân

Như chỉ ra trên hình trên, MatLab vẽ đường nội suy tuyến tính các điểm dữ liệu. Để đánh giá nhiệt độ tại điểm bất kỳ, chúng ta có thể dùng trực tiếp đồ thị. Ngược lại, ta có thể dùng hàm interp1:

>>t= interp1(hours, temps,9.3) % tính nhiệt độ tại giờ hours = 9.3
f=
22.9000
>> t= interp1(hours, temps,4.7) % tính nhiệt độ tại giờ hours =4.7
1=
22
>>t= interp1(hours, temps,[3.2 6.5 7.1 11.7]) % % tỉnh t° tại các giờ 3.2, 6.5...

t=
- 10.2000
30.0000
30.9000
24.9000

Đây là cách dùng ngầm định của hàm interp1. Cách viết interpl(x,y,x) có nghĩa là x- biến độc lập (hoành độ), y-biến phụ thuộc (tung độ), còn xạ là mảng chứa các giá trị nội suy. Ngoài ra còn ngầm định là nội suy tuyến tính.
Bây giờ thay cho đường thẳng nối các điểm dữ liệu ta giả thiết là dùng đường cong trơn nào đó để khớp các điểm số liệu. Giả thiết thông dụng nhất là dùng đa thức bậc 3, tức đa thức lập phương để mô phỏng từng đoạn giữa các điểm. số liệu liên tiếp sao cho độ dốc của từng đa thức bậc 3 trùng với điểm số liệu. Cách nội suy như vậy được gọi là chốt bậc 3 hay đơn giản nhất là chốt (spline). Chúng ta tìm lời giải này :

>>t= interp1(hours, temps,9.3,’spline’) % tính t tại giờ = 9.3
t=
21.8577
>> t= interp1(hours, temps,4.7, ’spline') % tính t" tại giờ =4,7
t=
22.3143
>>t= interp1(hours, temps,[3.2 6.5 7.1 11.7], 'spline')
t=
9.6734
30.0427
31.1755
25.3820

Hãy để ý là kết quả nội suy chốt có khác với nội suy tuyến tính. Do nội suy là quá trình đánh giá hoặc ước lượng giá trị nên ta cảm thấy nếu áp dụng các quy tắc ước lượng khác nhau sẽ dẫn đến các kết quả khác nhau.
Một cách dùng nội suy chốt thông dụng nhất là để làm mịn dữ liệu, có nghĩa là với một dãy số liệu cho trước ta dùng nội suy chốt để tính số liệu trong khoảng giữa thứ n. Ví dụ :

>> h= 1:0.1:12; % cứ 1 phần 10 tiếng lại đánh giá nhiệt độ một lần 
>> t = interp1(hours, temps,h, 'spline');
>> plot(hours, temps,'-', hours, temps,'+',h,t)
>> xlabel('hours'), ylabel(' Nhiet do Xen-so ')
>> title(' Nhiet do o Xuan Truong')

Trong hình trên, đường chấm chấm là đường nội suy tuyến tính, đường liền nét là nội suy chốt trơn phẳng, còn dữ liệu gốc được đánh dấu ‘+’. Bằng cách đòi hỏi độ phân giải mịn hơn trên trục giờ (hours) và dùng nội suy chốt ta có một đường cong mịn hơn, nhưng không phải là chính xác hơn, để tính nhiệt độ. Đặc biệt, độ dốc của lời giải chốt không thay đổi đột ngột tại các điểm dữ liệu.
Trước khi xét đến nội suy 2 chiều rất quan trọng là cần nhận ra hai hạn chế của interpl. Thứ nhất là không thể hỏi kết quả nằm ngoài phạm vi biến độc lập, tức interp1((hours, temps, 13.5) sẽ dẫn đến sai lỗi vì biến hours chỉ nằm trong khoảng giữa 1 và 12. Thứ hai là biến độc lập cần phải biến đổi đơn điệu, tức biến độc lập cần phải luôn luôn tăng hoặc luôn luôn giảm. Trong ví dụ của chúng ta hours là đơn điệu. Tuy nhiên nếu chúng ta xác định biến độc lập theo trình tự thời gian của ngày :

>>t_cua_ngay={7:12 1:6)%ngày bắt đầu lúc 7 giờ sáng và kết thúc lúc 6 giờ tối

t_cua_ngay =
7 8 9 10 11 12 1 2 3 4 5 6

thì biến độc lập không còn là đơn điệu nữa vì t_cua_ngay tăng tới 12 sau đó giảm xuống 1, sau đó lại tăng. Nếu dùng t_cua_ngay vào biến hours (giờ) trong interp1 thì sẽ xuất hiện lỗi. Cũng do nguyên nhân này nên không thể sử dụng temps để tìm giờ xảy ra nhiệt độ nào đó vì temps không đơn điệu.

3. Nội suy bậc 2

Mời các bạn đón xem phần 3.