Kiến thức

Vật lý tính toán – Phần 4: Nội suy đa thức – Mai trời sáng !

Vật lý tính toán – Phần 4: Nội suy đa thức

Nội suy đa thức (Polyminal Interpolation) là làm sao suy đoán ra hàm số mà đi qua gần đúng tất cả các điểm dữ liệu cho trước, nhờ đó chúng ta nhận xét được tính quy luật của chúng và dự đoán được điểm dữ liệu tiếp theo.

Đầu tiên chúng ta cùng nghiên cứu 2 phương pháp giúp nội suy đa thức có dạng thuần nhất y_n(x)= a_0+a_1x+a_2x^2+cdots+a_nx^n , ko chứa các hàm toán học đặc biệt như lượng giác hay logarit. n+1 là số điểm dữ liệu mà bạn có. Điều này nghe có vẻ kỳ lạ, ví dụ bạn khảo sát một đại lượng y mà bạn tin rằng nó quan hệ với đại lượng x theo hàm bậc 2, mà bạn đo n+1 lần, hổng nhẽ hàm của bạn phải là hàm bậc n? Thực ra hàm nội suy ko phản ánh đúng mối quan hệ thực tế của 2 đại lượng. Càng thu thập nhiều điểm dữ liệu thì hàm nội suy càng tiệm cận với hàm thực, các hệ số của các số hạng bậc cao sẽ dần bé đi thôi.

Nội suy Lagrange:

Đa thức nội suy có dạng: displaystyle y=sum^{n}_{i=0} L_i y_i trong đó displaystyle L_i= prod^{n}_{substack{j=0\jneq i}} frac{x-x_j}{x_i-x_j} .

 #include‹stdio.h› int main() { int n, i, j; printf("Nhap so diem cho truoc: "); scanf("%d",&n); double x[n], y[n], X, Y=0, L; printf("Nhap toa do cac diem:n"); for(i=0;i‹n;i++) { printf("(x%d,y%d)= ",i,i); scanf("%lf %lf",&x[i],&y[i]); } printf("X= "); scanf("%lf",&X); for(i=0;i‹n;i++) { L=1; for(j=0;j‹n;j++) if(j!=i) L*=(X-x[j])/(x[i]-x[j]); Y+=L*y[i]; } printf("Y= %lfn",Y); } 

Code:

computational-physics/Lagrange_Interpolation.c

Nhược điểm của phương pháp nội suy Lagrange là nếu bổ sung thêm điểm dữ liệu mới là phải tính lại từ đầu. Bây giờ là làm sao khi thêm dữ liệu mới vẫn tận dụng được đa thức đã tính trước đó. Chúng ta cùng bước sang phương pháp hiệu quả hơn, đó là nội suy Newton.

Nội suy Newton:

y= y_0+ dfrac{y_1-y_0}{x_1-x_0}(x-x_0) + dfrac{frac{y_2-y_1}{x_2-x_1}-frac{y_1-y_0}{x_1-x_0}}{x_2-x_0}(x-x_0)(x-x_1) + cdots
= y_0+ y'(x_0),(x-x_0) + y''(x_0),(x-x_0)(x-x_1) + cdots
trong đó y'(x_0), y''(x_0), ldots là các đạo hàm tại x=x0.
y'(x_0) = dfrac{y_1-y_0}{x_1-x_0} ,
y''(x_0) = dfrac{y'(x_1)-y'(x_0)}{x_2-x_0} , …
Chúng ta đặt ký hiệu như thế này: y= y[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)(x-x_1) + cdots
displaystyle Rightarrow y= sum^{n}_{i=0} f^{(i)}[x_0,x_1,ldots,x_i] prod^{i-1}_{j=0} (x-x_j) .
Tính các đạo hàm bằng cách tạo dựng ma trận D như sau:
D= begin{bmatrix} y_0 \ y_1 & f[x_0,x_1] \ y_2 & f[x_1,x_2] & f[x_0,x_1,x_2] \ vdots \ y_n & f[x_{n-1},x_n] & cdots & f[x_0,ldots,x_n] end{bmatrix}
Đạo hàm bậc cao được tạo nên bởi 2 đạo hàm bậc thấp hơn. Có cột đầu tiên là các giá trị đã biết, các cột tiếp theo được tính bằng cách D_{ij}=(D_{i,j-1}-D_{i-1,j-1})/(x_i-x_{i-j}) .

 #include‹stdio.h› int main() { int n, i, j; printf("Nhap so diem cho truoc: "); scanf("%d",&n); double x[n], y[n], X, Y=0, L, D[n][n]; for(i=0;i‹n;i++) { printf("(x%d,y%d)= ",i,i); scanf("%lf %lf",&x[i],&y[i]); } printf("X= "); scanf("%lf",&X); for(j=0;j‹n;j++) { if(j==0) for(i=0;i‹n;i++) D[i][j]=y[i]; else for(i=j;i‹n;i++) D[i][j]=(D[i][j-1]-D[i-1][j-1])/(x[i]-x[i-j]); L=1; for(i=0;i‹j;i++) L*=X-x[i]; Y+=D[j][j]*L; } printf("Y= %lfn",Y); } 

Code:

computational-physics/Newton_Interpolation.c

Cả 2 phương pháp ở trên đều có nhược điểm là số bậc của đa thức nội suy gắn liền với số lượng điểm dữ liệu thu thập. Tức là mặc dù càng nhiều thông tin thì việc dự đoán điểm tiếp theo càng chính xác, nhưng hàm nội suy sẽ càng phức tạp lên ko cần thiết, tai hại hơn là làm tăng khối lượng tính toán. Phương pháp nội suy Spline ra đời, giúp chúng ta tìm ra hàm nội suy với bậc thấp hơn và ko phụ thuộc vào số lượng điểm dữ liệu. Chi tiết về thuật toán này mình sẽ tìm hiểu sau.

Nội suy Lagrange, Newton, và cả nội suy Spline đều cố vẽ ra một hàm đi qua hết tất cả các điểm, điều này ko hay lắm vì trong các thí nghiệm vật lý, sai số đo đạc là điều tất yếu, các điểm dữ liệu đâu hoàn toàn nằm trên đường đồ thị lý thuyết mà quanh quẩn gần đường lý tưởng đó thôi. Càng cố nối điểm, chúng ta càng xa rời với đường lý tưởng hơn. Bởi vậy phương pháp bình phương tối thiểu được xem là hiệu quả nhất trong các thí nghiệm thực tế, đó là vẽ ra một đường sao cho tổng bình phương khoảng cách từng điểm dữ liệu tới đường đồ thị đó là nhỏ nhất, tức là thoả mãn điều kiện sum^{n}_{i=0} [y_i- sum^{m}_{j=0}a_j {x_i}^j]^2 ;mathrm{min} với m là số bậc của đa thức, ko phụ thuộc vào số điểm dữ liệu n.
displaystyle sum^{n}_{i=0} left[y_i- sum^{m}_{j=0}a_j {x_i}^jright]^2 ;mathrm{min}iffdisplaystyle frac{partial}{partial{a_k}} left[y_i- sum^{m}_{j=0}a_j {x_i}^jright]^2 = 0 với k = 0, 1, …, m .
iffdisplaystyle sum^{n}_{i=0} left[ -y_i{x_i}^k +sum^{m}_{k=0} a_j{x_i}^{j+k}right] = 0iffdisplaystyle sum^{n}_{i=0}sum^{m}_{k=0} a_j{x_i}^{j+k} = sum^{n}_{i=0} y_i{x_i}^k .
Ta biểu diễn được phương trình trên bằng phép toán ma trận như sau:
F^T Fa=F^T yiffbegin{bmatrix} 1 & 1 & 1 & cdots & 1 \ x_0 & x_1 & x_2 & cdots & x_n \ vdots \ x_0^m & x_1^m & x_2^m & cdots & x_n^m end{bmatrix} begin{bmatrix} 1 & x_0 & x_0^2 & cdots & x_0^m \ 1 & x_1 & x_1^2 & cdots & x_1^m \ vdots \ 1 & x_n & x_n^2 & cdots & x_n^m end{bmatrix} begin{bmatrix} a_0 \ a_1 \ vdots \ a_m end{bmatrix} = begin{bmatrix} 1 & 1 & 1 & cdots & 1 \ x_0 & x_1 & x_2 & cdots & x_n \ vdots \ x_0^m & x_1^m & x_2^m & cdots & x_n^m end{bmatrix} begin{bmatrix} y_0 \ y_1 \ vdots \ y_n end{bmatrix}
iffFa=yiffbegin{bmatrix} 1 & x_0 & x_0^2 & cdots & x_0^m \ 1 & x_1 & x_1^2 & cdots & x_1^m \ vdots \ 1 & x_n & x_n^2 & cdots & x_n^m end{bmatrix} begin{bmatrix} a_0 \ a_1 \ vdots \ a_m end{bmatrix} = begin{bmatrix} y_0 \ y_1 \ vdots \ y_n end{bmatrix}
Như vậy chúng ta chỉ cần dùng các phương pháp giải hệ phương trình như khử Gauss-Jordan chẳng hạn, giải ma trận:
left[ begin{array}{ccccc|c} 1 & x_0 & x_0^2 & cdots & x_0^m & y_0 \ 1 & x_1 & x_1^2 & cdots & x_1^m & y_1 \ vdots &&&&& vdots \ 1 & x_n & x_n^2 & cdots & x_n^m & y_n end{array} right]longrightarrowbegin{bmatrix} a_0 \ a_1 \ vdots \ a_m end{bmatrix}

 #include‹stdio.h› int main() { int m, n, i, j, k; printf("Nhap so bac cua da thuc: m= "); scanf("%d",&m); printf("Nhap so diem du lieu cho truoc: n= "); scanf("%d",&n); double x[n], y[n], X, Y=0, F[m+2][n], a[m+1], p; printf("Nhap gia tri cac diem cho truoc:n"); for(i=0;i‹n;i++) { printf("(x%d,y%d)= ",i,i); scanf("%lf %lf",&x[i],&y[i]); F[i][0]=1; for(j=1;j‹m+1;j++) F[i][j]=F[i][j-1]*x[i]; F[i][m+1]=y[i]; } for(k=0;k‹m+1||k‹n-1;k++) for(i=k+1;i‹n;i++) { p=F[i][k]/F[k][k]; for(j=k;j‹m+2;j++) F[i][j]-=p*F[k][j]; } for(k=k-1;k›0;k--) for(i=k-1;i›=0;i--) { p=F[i][k]/F[k][k]; for(j=k;j‹m+2;j++) F[i][j]-=p*F[k][j]; } for(k=0;k‹m+1||k‹n-1;k++) a[k]=F[k][m+1]/F[k][k]; printf("Nhap X= "); scanf("%lf",&X); for(j=m;j›=0;j--) Y=Y*X+a[j]; printf("Y= %lfn",Y); } 

Code:

computational-physics/Least_square_Method.c

Trả lời

Hủy trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

Gravatar

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com

Đăng xuất

 / 

Thay đổi

 )

Google photo

Bạn đang bình luận bằng tài khoản Google

Đăng xuất

 / 

Thay đổi

 )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter

Đăng xuất

 / 

Thay đổi

 )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook

Đăng xuất

 / 

Thay đổi

 )

Hủy bỏ

Connecting to %s

This site uses Akismet to reduce spam.

Learn how your comment data is processed

.

Chuyên mục: Kiến thức

Related Articles

Trả lời

Email của bạn sẽ không được hiển thị công khai. Các trường bắt buộc được đánh dấu *

Check Also
Close
Back to top button