Hàm đa thức - Chương 4: Giải hệ phương trình đại số tuyến tính

nclude 
#include 
#include 
#define spt 10
void main()
 {
	float a[spt][2*spt];
	float b[spt];
	int i,j,k,n,m,t;
	float max,c;
	char tl;
	clrscr();
	printf("Cho so phuong trinh n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran a :\n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	printf("\n");
	printf("Ma tran a ma ban da nhap");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	t=1;
	flushall();
	while (t)
	 {
	printf("Co sua ma tran a khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("Cho chi so cot can sua : ");
	scanf("%d",&j);
	printf("a[%d][%d] = ",i,j);
	scanf("%f",&a[i][j]);
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	printf("Ma tran a\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	printf("Cho cac phan tu cua ma tran b : \n");
	for (i=1;i<=n;i++)
	 {
	printf("b[%d] = ",i);
	scanf("%f",&b[i]);
	 }
	printf("\n");
	printf("Ma tran b ma ban da nhap\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("b[%d] = %15.5f\n",i,b[i]);
	printf("\n");
	t=1;
	flushall();
	while (t)
	 {
	printf("Co sua ma tran b khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("b[%d] = ",i);
	scanf("%f",&b[i]);
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	printf("\n");
	printf("Ma tran b\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("%15.5f\n",b[i]);
	printf("\n");
	t=1;
	flushall();
	i=1;
	while (t)
	 {
	if (a[i][i]==0)
	 {
	max=0;
	m=i;
	for (k=i+1;k<=n;k++)
	 if (max<fabs(a[k][i]))
	{
	 m=k;
	 max=fabs(a[i][i]);
	}
	 if (m!=i)
	{
	 for (j=i;j<=n;j++)
	{
	 c=a[i][j];
	 a[i][j]=a[m][j];
	 a[m][j]=c;
	}
	 c=b[i];
	 b[i]=b[m];
	 b[m]=c;
	}
	 if (m==i)
	{
	 t=0;
	 printf("MA TRAN SUY BIEN");
	}
	 }
	if (a[i][i]!=0)
	 {
	c=1/a[i][i];
	for (j=i;j<=n;j++)
	 a[i][j]=a[i][j]*c;
	b[i]=b[i]*c;
	for (k=1;k<=n;k++)
	 if (k!=i)
	{
	 c=a[k][i];
	 for (j=i;j<=n;j++)
	a[k][j]=a[k][j]-a[i][j]*c;
	 b[k]=b[k]-b[i]*c;
	}
	 }
	i=i+1;
	if (i==(n+1))
	 t=0;
	 }
	if (i==(n+1))
	 {
	printf("NGHIEM CUA HE");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("x[%d] = %15.5f\n",i,b[i]);
	 }
	getch();
 }
§3. PHƯƠNG PHÁP CHOLESKY
	Trong phương pháp Cholesky một ma trận đối xứng A được phân tích thành dạng A = RTR trong đó R là một ma trận tam giác trên.Hệ phương trình lúc đó chuyển thành AX = RTRX = B. Như vậy trước hết ta phân tích ma trận A thành tích hai ma trận. Sau đó giải hệ phương trình RTY = B và cuối cùng là hệ RX = Y. Chương trình mô tả thuật toán này được cho dưới đây:
Chương trình 4-5
#include 
#include 
#include 
#include 
#include 
#define max 6
void main()
 {
	float a[max][max],r[max][max];
	float b[max],x[max],y[max];
	int i,j,k,l,n,t;
	float s;
	char tl;
	clrscr();
	printf("Cho so phuong trinh n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran a : \n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	printf("\n");
	printf("Ma tran a ma ban da nhap\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	flushall();
	t=1;
	while (t)
	 {
	printf("Co sua ma tran a khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("Cho chi so cot can sua : ");
	scanf("%d",&j);
	printf("a[",i,",",j,"] = ");
	scanf("%f",&a[i][j]);
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	printf("Ma tran a\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	printf("Cho cac phan tu cua ma tran b : \n");
	for (i=1;i<=n;i++)
	 {
	printf("b[%d] = ",i);
	scanf("%f",&b[i]);
	 }
	printf("\n");
	printf("Ma tran b ma ban da nhap\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("b[%d] = %15.5f\n",i,b[i]);
	printf("\n");
	flushall();
	t=1;
	while (t)
	 {
	printf("Co sua ma tran b khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("b[%d] = ",i);
	scanf("%f",&b[i]);
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	printf("\n");
	printf("Ma tran b\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("b[%d] = %15.5f\n",i,b[i]);
	printf("\n");
	for(i=1;i<=n;i++)
	 for(j=1;j<=n;j++)
	r[i][j]=0.0;
	for (i=1;i<=n;i++)
	 {
	if (a[i][i]>=0)
	 {
	r[i][i]=sqrt(a[i][i]);
	for (j=1+i;j<=n;j++)
	 r[i][j]=a[i][j]/r[i][i];
	for (k=i+1;k<=n;k++)
	 for (l=k;l<=n;l++)
	a[k][l]=a[k][l]-r[i][k]*r[i][l];
	 }
	 }
	for (k=1;k<=n;k++)
	 {
	s=b[k];
	if (k!=1)
	 for (i=1;i<=k-1;i++)
	s=s+r[i][k]*y[i];
	 y[k]=-s/r[k][k];
	 }
	for (i=n;i>=1;i--)
	 {
	s=-y[i];
	if (i!=n)
	 for (k=i+1;k<=n;k++)
	s=s-r[i][k]*x[k];
	 x[i]=s/r[i][i];
	 }
	printf("Nghiem cua he phuong trinh la\n ");
	for (i=1;i<=n;i++)
	 printf("x[%d] = %10.5f\n",i,x[i]);
	getch();
 }
§4. PHƯƠNG PHÁP CROUT
	Phương pháp Crout là một dạng của phương pháp Gauss.Với phương pháp Gauss, chúng ta biến đổi ma trận A thành một ma trận tam giác thì ở phương pháp Crout chúng ta phân tích ma trận này thành tích của ma trận tam giác trên R và ma trận tam giác dưới L. Trong ma trận L, các hệ số trên đường chéo chính bằng 1. Như vậy phương trình AX = B được viết thành :
	A.X = L.R.X = B
Chúng ta đặt:
	RX = Y
nên :	LY = B
Như vậy trước hết chúng ta phân tích ma trận thành tích của L.R. Tiếp theo ta giải phương trình LY = B và sau đó giải phương trình RX = A để tìm nghiệm X.
Chương trình 4-6
#include 
#include 
#include 
#include 
#include 
#define max 6
void main()
 {
	float b[max],x[max],y[max];
	float a[max][max],r[max][max],l[max][max];
	int i,j,k,n,t;
	float c,tr,tl,s;
	char tloi;
	clrscr();
	printf("Cho so phuong trinh n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran a : \n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	printf("\n");
	printf("Ma tran a ma ban da nhap");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	printf("%10.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	t=1;
	flushall();
	while (t)
	 {
	printf("Co sua ma tran a khong(c/k)?");
	scanf("%c",&tloi);
	if (toupper(tloi)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("Cho chi so cot can sua : ");
	scanf("%d",&j);
	printf("a[%d][%d] = ",i,j);
	scanf("%f",&a[i][j]);
	flushall();
	 }
	if (toupper(tloi)=='K')
	 t=0;
	 }
	printf("Ma tran a\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%10.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	printf("Cho cac phan tu cua ma tran b : \n");
	for (i=1;i<=n;i++)
	 {
	printf("b[%d] = ",i);
	scanf("%f",&b[i]);
	 }
	printf("\n");
	printf("Ma tran b ma ban da nhap");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("b[%d] = %10.5f\n",i,b[i]);
	printf("\n");
	t=1;
	flushall();
	while (t)
	 {
	printf("Co sua ma tran b khong(c/k)?");
	scanf("%c",&tloi);
	if (toupper(tloi)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("b[%d] = ",i);
	scanf("%f",&b[i]);
	flushall();
	 }
	if (toupper(tloi)=='K')
	 t=0;
	 }
	printf("\n");
	printf("Ma tran b\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("b[%d] = %10.5f\n",i,b[i]);
	printf("\n");
	for(i=1;i<=n;i++)
	 for(j=1;j<=n;j++)
	{
	 r[i][j]=0.0;
	 l[i][j]=0.0;
	}
	for (i=1;i<=n;i++)
	 {
	r[1][i]=a[1][i];
	l[i][i]=1.0;
	l[i][1]=a[i][1]/a[1][1];
	 }
	for (k=2;k<=n;k++)
	 {
	for (j=k;j<=n;j++)
	 {
	tr=0.0;
	for (i=1;i<=k;i++)
	 tr=tr+l[k][i]*r[i][j];
	r[k][j]=a[k][j]-tr;
	 }
	if (k!=n)
	 {
	for (i=1;i<=n;i++)
	 {
	tl=0.0;
	for (j=1;j<=k-1;j++)
	 tl=tl+l[i][j]*r[j][k];
	l[i][k]=(a[i][k]-tl)/r[k][k];
	 }
	 }
	else
	 printf("\n");
	 }
	if (l[1][1]==0.0)
	 if (b[1]==0.0)
	printf("He da cho vo nghiem\n");
	 else
	{
	 printf("He da cho co vo so nghiem\n");
	 y[n]=c;
	}
	else
	 y[1]=b[1]/l[1][1];
	for (i=2;i<=n;i++)
	 {
	s=0.0;
	for (k=1;k<=i-1;k++)
	 s=s+l[i][k]*y[k];
	y[i]=(b[i]-s)/l[i][i];
	 }
	if (r[n][n]==0.0)
	 if (y[n]==0.0)
	printf("He da cho vo nghiem\n");
	 else
	{
	 printf("He da cho co vo so nghiem\n");
	 x[n]=c;
	}
	else
	 x[n]=y[n]/r[n][n];
	for (i=n-1;i>=1;i--)
	 {
	s=0.0;
	for (k=i+1;k<=n;k++)
	 s+=r[i][k]*x[k];
	x[i]=(y[i]-s)/r[i][i];
	 }
	printf("\n");
	printf("Nghiem cua he da cho la\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 printf("x[%d] = %15.5f\n",i,x[i]);
	getch();
 }
§5. PHƯƠNG PHÁP LẶP ĐƠN
Xét hệ phương trình AX = F. Bằng cách nào đó ta đưa hệ phương trình về dạng 
	X = BX + G
trong đó:
	B = (bij)n,n
	G = (g1,g2,...,gn)T
Chọn vectơ:
	X = ( x1(o),x2(o),....,xn(o) )T 
làm xấp xỉ thứ 0 của nghiệm đúng và xây dựng xấp xỉ 
	X(m+1) = BX(m) + G ( m = 0,1,....)
Người ta chứng minh rằng nếu phương trình ban đầu có nghiệm duy nhất và một trong ba chuẩn của ma trận B nhỏ hơn 1 thì dãy xấp xỉ hội tụ về nghiệm duy nhất đó.(Cho một ma trận B,chuẩn của ma trận B, kí hiệu là || B || là một trong 3 số :
(Chuẩn của ma trận quan hệ tới sự hội tụ của phương pháp lặp)
Ví dụ: Chúng ta có phương trình
Chúng ta đưa phương trình về dạng :
Như vậy : 
	 và 
Dễ thấy ; và nên phép lặp hội tụ. Chương trình lặp đơn là:
Chương trình 4-7
#include 
#include 
#include 
#include 
#include 
#define max 10
void main()
 {
	float a[max][max];
	float f[max],x0[max],x1[max];