Giải hệ phương trình: Code: computational-physics/Gauss_Elimination_with_partial_pivoting.c Code: computational-physics/Gauss-Jordan_Elimination.c
Mình lược bỏ Interchange cho dễ nhìn. Bạn có thể thấy code khử Gauss và khử Jordan tương đồng nhau, chẳng qua chiều con chạy i, j chạy ngược nhau.
Sở dĩ có bước tráo hàng [interchange] dòng 14-20 trong bước khử Gauss là để tối ưu thuật toán, bởi phép chia cho một số lớn thì ít sai số hơn so với chia cho một số bé, đặc biệt là tránh chia cho số cực bé. Với bài toán đơn giản, ko quá cần thiết phải có bước tráo hàng này. Code ở trên chỉ thực hiện tráo một phần, nếu cầu toàn hơn, bạn có thể tráo toàn phần, tức là vừa tráo hàng vừa tráo cột. Bạn nào yêu thích lập trình có thể thử, chứ cũng ko cần thiết phải biết.
Cách 1:
#include‹stdio.h›
int main[] {
int n, i, j, k;
printf["Nhap co ma tran cua Ax=b\nn= "]; scanf["%d",&n];
double a[n][n+1], p;
printf["Nhap gia tri cho ma tran:\n"];
for[i=0;i‹n;i++] {
printf["Hang thu %d: ",i+1];
for[j=0;j‹n+1;j++] scanf["%lf",&a[i][j]];
}
/*Gauss Elimination*/
for[k=0;k‹n-1;k++] for[i=k+1;i‹n;i++] {
p=a[i][k]/a[k][k];
for[j=k;j‹n+1;j++] a[i][j]-=p*a[k][j];
}
/*Jordan Elimination*/
for[k=n-1;k›0;k--] for[i=k-1;i›=0;i--] {
p=a[i][k]/a[k][k];
for[j=k;j‹n+1;j++] a[i][j]-=p*a[k][j];
}
/*Show solution*/
printf["Nghiem he phuong trinh la:\n"];
for[i=0;i‹n;i++] printf["x%d= %.3lf\n",i+1,a[i][n]/a[i][i]];
}
Cách 2:
Code: computational-physics/Synchronic_Gauss-Jordan_Elimination.c
Chỉ sai khác một chút xíu mà nghiệm của 2 hệ hoàn toàn lệch nhau [bởi sau khi khử hàng, nó để lại số hạng cực bé, như mình trình bày vì sao có chèn thêm bước interchange trong thuật toán]. Nếu như chuyện sai số từ dữ liệu đầu vào là ko thể tránh khỏi, vậy thì hãy chấp nhận sai số trong giới hạn bao nhiêu chữ số thập phân. Phương pháp lặp là phương pháp thử giá trị cho nghiệm, qua các bước lặp, các nghiệm giả định sẽ dần hội tụ về nghiệm chuẩn xác của hệ [dĩ nhiên trong phạm vi sai số bạn chọn].
Tuy nhiên, phương pháp lặp chỉ hữu hiệu hơn phương pháp khử hàng khi xử lý các hệ có tính ổn định kém, chứ thực tình trong nhiều bài toán, phương pháp lặp còn “phế” hơn nhiều, bởi nó vướng phải một điểm yếu, đó là để đảm bảo sự hội tụ của nghiệm, hệ phải thoả mãn điều kiện tính chéo trội, tức là độ lớn của phần tử trên đường chéo chính phải lớn hơn độ lớn của các phần tử trong cùng hàng cộng lại.
#include‹stdio.h› int main[] { int n, i, j, k=0; printf["Nhap co he phuong trinh:\nn= "]; scanf["%d",&n]; double a[n][n+1]; printf["Nhap gia tri cho ma tran:\n"]; for[i=0;i‹n;i++] { printf["Hang thu %d: ",i+1]; for[j=0;j‹n+1;j++] scanf["%lf",&a[i][j]]; } double x0[n], x[n], S, err=1e-9, delta; for[i=0;i‹n;i++] x0[i]=0; do { for[i=0;i‹n;i++] { S=0; for[j=0;j‹n;j++] if[j!=i] S+=a[i][j]*x0[j]; x[i]=[a[i][n]-S]/a[i][i]; } delta=0; for[i=0;i‹n;i++] { delta+=[x[i]-x0[i]]*[x[i]-x0[i]]; x0[i]=x[i]; } k++; } while[delta›err]; printf["Nghiem he pt la:\n"]; for[i=0;i‹n;i++] printf["x[%d]= %lf\n",i+1,x[i]]; printf["Voi so buoc lap la k= %d\n",k]; }Code: computational-physics/Jacobi_Iteration.c
Dòng 12 giả định tập nghiệm có giá trị bằng 0 hết [đặt bằng bao nhiêu cũng dc]. Dòng 14-18 chính là bước lặp theo công thức . Dòng 19-23 là tính sai số delta bằng tổng bình phương của các hiệu giữa nghiệm bước lặp cũ và nghiệm bước lặp mới, sao cho đến khi nhỏ hơn err cho trước thì dừng [1e-9 nghĩa là 1*10^{-9}].
Code: computational-physics/Gauss-Seidel_Iteration.c
Tìm ma trận nghịch đảo:
Áp dụng phương pháp khử Gauss-Jordan để biến đổi ma trận dạng [A|I] thành dạng [I|B], khi đó B là ma trận nghịch đảo của A.
Code: computational-physics/Inverse_Matrix.c
Dòng 13-16 để kiểm tra phần tử a[k][k], nếu nó bằng 0 thì khẳng định ma trận A ko thể nghịch đảo luôn, vì bắt buộc để nửa trái trở thành ma trận đơn vị, phần tử đó phải là pivot. Lưu ý rằng, ta chỉ được phép khử hàng mà ko được phép tráo hàng, vì như vậy sẽ thành một ma trận khác ko tương đương, hay chính xác thì tráo hàng chỉ được dùng khi giải hệ phương trình [bởi khi tráo hàng ta tráo luôn cả phần tử b của phương trình Ax=b, điều đó ko làm ảnh hưởng đến vector nghiệm x, nhưng thực chất ma trận A đã bị biến đổi].
Tính định thức:
Như mình đã nói ở bài viết trong seri Đại số tuyến tính, ko ai đi dùng định nghĩa để tính định thức cả, mà biến đổi ma trận A thành ma trận tam giác trên U [mối quan hệ giữa định thức của 2 ma trận tuân theo tính chất 1, 2, 3 cũng đều nêu trong bài viết đó], và định thức của U bằng tích các phần tử đường chéo chính.
#include‹stdio.h› #include‹math.h› int main[] { int n, i, j, k, imax, sign=1; printf["Nhap co ma tran vuong n*n:\nn= "]; scanf["%d",&n]; double a[n][n], p; printf["Nhap gia tri cho ma tran:\n"]; for[i=0;i‹n;i++] { printf["Hang thu %d: ",i+1]; for[j=0;j‹n;j++] scanf["%lf",&a[i][j]]; } for[k=0;k‹n-1;k++] { imax=k; for[i=k+1;i‹n;i++] if[fabs[a[i][k]]›fabs[a[imax][k]]] imax=i; if[imax!=k] sign*=-1; for[j=0;j‹n;j++] { p=a[k][j]; a[k][j]=a[imax][j]; a[imax][j]=p; } for[i=k+1;i‹n;i++] { p=a[i][k]/a[k][k]; for[j=k;j‹n;j++] a[i][j]-=p*a[k][j]; } } p=1; for[i=0;i‹n;i++] p*=a[i][i]; printf["Det= %.3lf\n",sign*p]; }Code: computational-physics/Determinant.c
Sở dĩ dòng 15 có biến sign vì theo tính chất 2 của định thức, tráo hàng ma trận là định thức của ma trận đó đổi dấu, nên sign sẽ lưu dấu của định thức mỗi khi xảy ra tráo hàng, cuối cùng mới nhân vào tích.