Hàm đa thức - Chương 2: Giải gần đúng phương trình đại số và siêu việt

("Cho khoang can tim nghiem [a,b]\n");
	printf("Cho gia tri duoi a = ");
	scanf("%f",&x2);
	printf("Cho gia tri tren b = ");
	scanf("%f",&x1);
	if (f(x1)*f(x2)>0)
	 {
	printf("\n");
	printf("Nghiem khong nam trong doan nay\n");
	getch();
	exit(1);
	 }
	eps=1e-5;
	x0=(x1+x2)/2;
	dem=0;
	do
	 {
	dem=dem+1;
	h1=x1-x0;
	h2=x0-x2;
	gamma=h2/h1;
	a=(gamma*f(x1)- f(x0)*(1+gamma)+f(x2))/(gamma*(h1*h1)*(1+gamma));
	b=(f(x1)-f(x0)-a*(h1*h1))/h1;
	c=f(x0);
	if ((a==0)&&(b!=0))
	 {
	n1=-c/b;
	n2=n1;
	 }
	if ((a!=0)&&(b==0))
	 {
	n1=(-sqrt(-c/a));
	n2=(sqrt(-c/a));
	 }
	if ((a!=0)&&(b!=0))
	 {
	n1=x0-2*c/(b+(sqrt(b*b-4*a*c)));
	n2=x0-2*c/(b-(sqrt(b*b-4*a*c)));
	 }
	if (fabs(n1-x0)>fabs(n2-x0))
	 xr=n2;
	else
	 xr=n1;
	if (xr>x0)
	 {
	x2=x0;
	x0=xr;
	 }
	else
	 {
	x1=x0;
	x0=xr;
	 }
	 }
	while (fabs(f(xr))>=eps);
	printf("\n");
	printf("Phuong trinh co nghiem x = %.5f sau %d lan lap",xr,dem);
	getch();
 }
float f(float x)
 {
	float a=sin(x)-x/2;
	return(a);
 }
§7. PHƯƠNG PHÁP LẶP BERNOULLI
	Có nhiều phương pháp để tìm nghiệm của một đa thức. Ta xét phương trình:
 	aoxn + a1xn-1 + ××× + an = 0
Nghiệm của phương trình trên thoả mãn định lí: Nếu max{| a1 |, | a2 |,..., |an |} = A thì các nghiệm của phương trình thoả mãn điều kiện | x | < 1 + A/ | a0|
 	Phương pháp Bernoulli cho phép tính toán nghiệm lớn nhất a của một đa thức Pn(x) có n nghiệm thực phân biệt. Sau khi tìm được nghiệm lớn nhất a ta chia đa thức Pn(x) cho (x-a) và nhận được đa thức mới Qn-1(x). Tiếp tục dùng phương pháp Bernoulli để tìm nghiệm lớn nhất của Qn-1(x).
 Sau đó lại tiếp tục các bước trên cho đến khi tìm hết các nghiệm của Pn(x).
	Chúng ta khảo sát phương trình sai phân j có dạng như sau :
	j = aoyk+n + a1yk+n-1 +.....+ anyk = 0	(1)
Đây là một phương trình sai phân tuyến tính hệ số hằng. Khi cho trước các giá trị đầu yo, y1,..yn-1 ta tìm được các giá trị yn, yn+1,.. Chúng được gọi là nghiệm của phương trình sai phân tuyến tính (1).
Đa thức 
Pn(x) = a0xn + a1xn-1 +..+an-1x + an 	(2)
với cùng một hệ số ai như (1) được gọi là đa thức đặc tính của phương trình sai phân tuyến tính (1). Nếu (2) có n nghiệm phân biệt x1, x2,.., xn thì (1) có các nghiệm riêng là 
Nếu yi là các nghiệm của phương trình sai phân là tuyến tính (1),thì 
	(3)
với các hệ số ci bất kì cũng là nghiệm của phương trình sai phân tuyến tính hệ số hằng (1).
Nếu các nghiệm cần sao cho :
	| x1| ³ | x2 | ³...| xn|
Vậy 	
và	
do đó : 
do 	x1 > x2 >...> xn
nên:	 khi k ® ¥ 
vậy:	 khi k ® ¥ 
Nghĩa là :	
	Nếu phương trình vi phân gồm n+1 hệ số, một nghiệm riêng yk có thể được xác định từ n giá trị yk-1, yk-2,...,yn-1. Điều cho phép tính toán bằng cách 
truy hồi các nghiệm riêng của phương trình vi phân.
	Để tính nghiệm lớn nhất của đa thức, ta xuất phát từ các nghiệm riêng y1 = 0, y1 = 0,.., yn =1 để tính yn+1. Cách tính này được tiếp tục để tính yn+2 xuất phát từ y1 = 0, y2 = 0,..,yn+1 và tiếp tục cho đến khi yk+1/yk không biến đổi nữa. Trị số của yk+n được tính theo công thức truy hồi :
	(4)
Ví dụ: Tính nghiệm của đa thức Pn(x) = P3(x) = x3 - 10x2 + 31x - 30. 
Như vậy ao = 1, a1 = -10,a2 = 31 và a3 = -30.
Phương trình sai phân tương ứng là :
yk+3 -10yk+2 + 31yk+1 - 30yk = 0
Ta cho trước các giá trị y1 = 0; y2 = 0 và y3 = 1. Theo (4) ta tính được :
y4 = - (-10y3 + 31y2 - 30y1) = 10
y5 = - (-10y4 + 31y3 - 30y2) = 69
y6 = - (-10y5 + 31y5 - 30y3) = 410
y7 = - (-10y6 + 31y5 - 30y4) = 2261
y8 = - (-10y7 + 31y6 - 30y5) = 11970
y9 = - (-10y8 + 31y7 - 30y6) = 61909
y10 = - (-10y9 + 31y8 - 30y8) = 315850
y11 = - (-10y10 + 31y9 - 30y8) = 1598421
y12 = - (-10y11 + 31y10 - 30y9) = 8050130
y13 = - (-10y12 + 31y11 - 30y10) = 40425749
y14 = - (-10y13 + 31y12 - 30y11) = 202656090
y15 = - (-10y14 + 31y13 - 30y12) = 1014866581
y16 = - (-10y15 + 31y14 - 30y13) = 5079099490
y17 = - (-10y16 + 31y15 - 30y14) = 24409813589
y18 = - (-10y17 + 31y16 - 30y15) = 127092049130
y19 = - (-10y18 + 31y17 - 30y16) = 635589254740
Tỉ số các số yk+1/yk lập thành dãy :
	10 ; 6.9 ; 5.942 ; 5.5146 ; 5.2941 ; 5.172 ; 5.1018 ; 5.0607 ; 5.0363 ; 5.0218 ; 5.013 ; 5.0078 ; 5.0047 ; 5.0028 ; 5.0017 ; 5.001 
nghĩa là chúng sẽ hội tụ tới nghiệm lớn nhất là 5 của đa thức.
Chương trình 2-7
//phuong phap Bernoulli
#include 
#include 
#include 
#include 
#define max 50
void main()
 {
	float a[max],y[max];
	int k,j,i,n,l;
	float s,e1,e2,x0,x1,x;
	clrscr();
	printf("Cho bac cua da thuc can tim nghiem n = ");
	scanf("%d",&n);
	e1=1e-5;
	printf("Cho cac he so cua da thuc can tim nghiem\n");
	for (i=0;i<=n;i++)
	 {
	printf("a[%d] = ",i);
	scanf("%f",&a[i]);
	 }
	for (k=0;k<=n;k++)
	 a[k]=a[k]/a[0];
	 tt: x1=0;
	for (k=2;k<=n;k++)
	 y[k]=0;
	y[1]=1;
	l=0;
	do
	 {
	l=l+1;
	s=0;
	for (k=1;k<=n;k++)
	 s=s+y[k]*a[k];
	y[0]=-s;
	x=y[0]/y[1];
	e2=fabs(x1 - x);
	x1=x;
	for (k=n;k>=1;k--)
	y[k]=y[k-1];
	 }
	while((l=e1));
	if(e2>=e1)
	 {
	printf("Khong hoi tu");
	getch();
	exit(1);
	 }
	else
	 printf("Nghiem x = %.4f\n",x);
	n=n-1;
	if (n!=0)
	 {
	a[1]=a[1]+x;
	for (k=2;k<=n;k++)
	 a[k]=a[k]+x*a[k-1];
	goto tt;
	 }
	getch();
 }
Kết quả nghiệm của đa thức x3 - 10x2 + 31x - 30 là:5 ; 3 và 2
§8. PHƯƠNG PHÁP LẶP BIRGE - VIETTE 
	Các nghiệm thực, đơn giản của một đa thức Pn(x) được tính toán khi sử dụng phương pháp Newton	
	(1)
Để bắt đầu tính toán cần chọn một giá trị ban đầu xo. Chúng ta có thể chọn một giá trị xo nào đó, ví dụ :
và tính tiếp các giá trị sau :	
Tiếp theo có thể đánh giá Pn(xi) theo thuật toán Horner :
P0 = a0
	P1 = P0xi + a1	 	(2)
	P2 = P1xi + a2
	P3 = P2xi + a3
 	..................
	P(xi) = Pn = Pn-1xi + an
Mặt khác khi chia đa thức Pn(x) cho một nhị thức (x - xi) ta được :
Pn(x) = (x - xi)Pn-1(x) + bn	(3)
với bn = Pn(xi). Đa thức Pn-1(x) có dạng:
	Pn-1(x) = boxn-1 + b1xn-2 + p3xn-3 +..+ bn-2x + bn-1	(4)
Để xác định các hệ số của đa thức (4) ta thay (4) vào (3) và cân bằng các hệ số với đa thức cần tìm nghiệm Pn(x) mà các hệ số ai đã cho:
	(x - xi)( boxn-1 + b1xn-2+b3xn-3 +..+ bn-2x + bn-1 ) + bn
	= aoxn + a1xn-1 + a2xn-2 +...+ an-1x + a=	(5)
Từ (5) rút ra :
bo = ao
b1 = a1 + boxi	(6)
	b2 = a2 + b1xi
	......
	bk = ak + bk-1xi
	.......
	bn = an + bn-1xi = Pn(xi)
	Đạo hàm (3) ta được :
và:	(7)
	Như vậy với một giá trị xi nào đó theo (2) ta tính được Pn(xi) và kết hợp (6) với (7) tính được P¢n(xi). Thay các kết quả này vào (1) ta tính được giá trị xi+1. Quá trình được tiếp tục cho đến khi | xi+1 - xi | < e hay Pn(xi+1) » 0 nên a1 » xi+1 là một nghiệm của đa thức.
	Phép chia Pn(x) cho (x - a1) cho ta Pn-1(x) và một nghiệm mới khác được tìm theo cách trên khi chọn một giá trị xo mới hay chọn chính xo=a1. Khi bậc của đa thức giảm xuống còn bằng 2 ta dùng các công thức tìm nghiệm của tam thức để tìm các nghiệm còn lại.
Ví dụ: Tìm nghiệm của đa thức P3(x) = x3 - x2 -16x + 24 
ao = 1	a1 = -1	a2= -16	a3 = 24
Chọn xo = 3.5 ta có :
Po = ao = 1
P1 = a1 + pox0 = -1 + 3.5*1 = 2.5	
P2 = a2 + p1x0 = -16 + 3.5*2.5 = -7.25
P3 = a3 + p2x0 = 24 + 3.5*(-7.25) = - 1.375
b0 = a0 = 1;
b1 = a1 + box0 = -1 + 3.5*1 = 2.5	
b2 = a2 + b1x0 = -16 + 3.5*2.5 = -7.25
P2(3.5) = b0x2 + b1x + b2 = 13.75
Lặp lại bước tính trên cho x1 ta có:
Po = ao = 1
P1 = a1 + pox1 = -1 + 3.6*1 = 2.6	
P2 = a2 + p1x1 = -16 + 3.6*2.6 = -6.64
P3 = a3 + p2x1 = 24 + 3.6*(-6.64) = - 0.096
bo = ao = 1
b1 = a1 + box1 = -1 + 3.6*1 = 2.6	
b2 = a2 + p1x1 = -16 + 3.6*2.6 = -6.64
P2(3.6) = b0x2 + b1x + b2 = 15.68
Quá trình cứ thế tiếp tục cho đến khi sai số chấp nhận được. Chương trình dưới đây mô tả thuật tính trên.
Chương trình 2-8
//phuong phap Birge-Viette
#include 
#include 
#include 
#define max 20
void main()
 {
	float a[max],p[max],d[max],x[max];
	int k,j,i,n;
	float e1,e2,x0,x1;
	clrscr();
	printf("Cho bac cua da thuc n = ");
	scanf("%d",&n);
	e1=0.0001;
	printf("Cho cac he so cua da thuc can tim nghiem\n");
	for (i=0;i<=n;i++)
	 {
	printf("a[%d] = ",i);
	scanf("%f",&a[i]);
	 }
	x0=a[0];
	for (i=0;i<=n;i++)
	 a[i]=a[i]/x0;
	printf("Nghiem cua phuong trinh : \n");
	tt:x0=-a[n]/a[n-1];
	j=0;
	do
	 {
	j=j+1;
	p[1]=x0+a[1];
	d[1]=1.0;
	for (k=2;k<=n;k++)
	 {
	p[k]=p[k-1]*x0+a[k];
	d[k]=d[k-1]*x0+p[k-1];
	 }
	x1=x0-p[n]/d[n];
	e2=fabs(x1-x0);
	if (e2>e1)
	x0=x1;
	 }
	while((j=e1));
if (e2>=e1)
	 printf("Khong hoi tu");
	else
	 printf(" x = %.4f\n",x1);
	n=n-1;
	if (n!=0)
	 {
	for (k=1;k<=n;k++)
	 a[k]=p[k];
	goto tt;
	 }
	getch();
 }
	Dùng chương trình trên để tìm nghiệm của đa thức x4 + 2x3 - 13x2 - 14x + 24 ta được các nghiệm là:-4 ; 3 ; -2 và 1.
§9. PHƯƠNG PHÁP NGOẠI SUY AITKEN
	Xét phương pháp lặp :
	x = f(x)	(1)
với f(x) thoả mãn điều kiện hội tụ của phép lặp, nghĩa là với mọi xÎ [a, b] ta có:
| f’(x) | £ q < 1	(2)
Như vậy :
	xn+1 = f(xn)	(3)
	xn = f(xn-1)	(4)
Trừ (3) cho (4) và áp dụng định lí Lagrange cho vế phải với c Î [a, b] ta có :
	xn+1- xn = f(xn) - f(xn-1) = (xn - xn-1)f’(c)	(5)
Vì phép lặp (1) nên :
	| xn+1- xn | £ q | xn - xn-1 |	(6)
b = an - pb-2
Chúng ta nhận thấy rằng a được tính toán xuất phát từ cùng một công thức truy hồi như các hệ số bk và tương ứng với hệ số bn-1
	bn-1 = an-1 + sbn-2 - pbn-3 = a
Hệ số bn là :
	bn = an + sbn-1 - pbn-2 = sbn-1 + b
và cuối cùng :
	R1(x) = ax + b = b+-1(x - s) + bn
	Ngoài ra các hệ số bi phụ thuộc vào s và p và bây giờ chúng ta cần phải tìm các giá trị đặc biệt s* và p* để cho bn-1 và bn triệt tiêu. Khi đó r1(x)= 0 và nghiệm của tam thức x2 - s*x + p*x sẽ là nghiệm của đa thức Pn(x). Ta biết rằng bn-1 và bn là hàm của s và p :
	bn-1 = f(s, p)
	bn = g(s, p)
Việc tìm s* và p* đưa đến việc giải hệ phương trình phi tuyến:
Phương trình này có thể giải dễ dàng nhờ phương pháp Newton. Thật vậy với một phương trình phi tuyến ta có công thức lặp:	
xi+1 = xi - f(xi)/f'(xi) 
hay	f'(xi)(xi+1 - xi) = -f(xi)
Với một hệ có hai phương trình,công thức lặp trở thành:
	J(Xi)(Xi+1 - Xi) = -F(Xi)