Giáo trình C - Chương 8: Giải gần đúng phương trình đại số và siêu việt

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) 
với Xi = { si,pi}
T và Xi+1 = { si+1,pi+1}
T
108
)p,s(g
)p,s(f
)X(F
ii
ii
i = 
iJ X
f
s
f
p
g
s
g
p
( )=
∂
∂
∂
∂
∂
∂
∂
∂
Quan hệ : J(Xi)∆X = -F(Xi) với ∆X = {si+1 - si,pi+1 - pi}T t−ơng ứng với một hệ ph−ơng trình 
tuyến tính hai ẩn số ∆s = si+1 - si và ∆p = pi+1 - pi : 
∂
∂
∂
∂
∂
∂
∂
∂
f
s s
f
p
p f s p
g
s s
g
p
p g s p
i i
i i
∆ ∆
∆ ∆
+ =
+ =
⎧
⎨
⎪⎪⎪
⎩
⎪⎪⎪
−
−
( )
( )
,
,
Theo công thức Cramer ta có : 
 ∆s
f
g
p
g
f
p=
− +∂∂
∂
∂
δ 
 ∆p
g
f
s
f
g
s=
− +∂∂
∂
∂
δ 
 δ ∂∂
∂
∂
∂
∂
∂
∂= −
f
s
g
p
f
p
g
s
Để dùng đ−ợc công thức này ta cần tính đ−ợc các đạo hàm ∂∂
f
s
,
∂
∂
f
p
,
∂
∂
g
s
,
∂
∂
g
p
.Các đạo hàm này 
đ−ợc tính theo công thức truy hồi. 
Do bo = ao nên 
 ob
s
∂
∂ =0
ob
p
∂
∂ =0
b1 = a1 + sbo nên 1∂
∂
b
s bo
= 1 0∂∂
b
p
= 
b2 = a2 + sb1- pbo nên 2 2 1∂∂
∂
∂
∂
∂
∂
∂
b
s
a
s
sb
s
pb
s
o= + −
( ) ( ) 
Mặt khác : 2 0∂∂
a
s
= 1 1 1
∂
∂
∂
∂
( )sb
s s
b
s b= +
 opb
s
∂
∂
( )
=0 
nên : 
 2
1
∂
∂
b
s b sbo
= + 
b3 = a3 + sb2- pb1 nên 
 3
2
2 1∂∂
∂
∂
∂
∂
b
s b s
b
s
p b
s
= + − 
Nếu chúng ta đặt : 
 k
k
b
s c
∂
∂ = −1
thì : co = bo (2) 
 c1 = b1 + sbo = b1 + sco 
109
 c2 = b2 + sc1 - pco 
 .................... 
 ck = bk + sck-1 - pck-2 
 cn-1 = bn-1 + scn-2 - pcn-3 
Nh− vậy các hệ số cũng đ−ợc tính theo cách nh− các hệ số bk.Cuối cùng với f = bn-1 và g = bn 
ta đ−ợc: 
2n1n3n2n cs
f
c
s
f
c
s
f
c
s
f
−−−− =∂
∂=∂
∂=∂
∂=∂
∂
2
2n3n1n
3nn2n1n
ccc
cbcb
s
−−−
−−−
−
−=∆ (3) 
2
2n3n1n
2nn1n1n
ccc
cbcb
p
−−−
−−−
−
−=∆ (4) 
 Sau khi phân tích xong Pn(x) ta tiếp tục phân tích Pn-2(x) theo ph−ơng pháp trên Các 
b−ớc tính toán gồm : 
 - Chọn các giá trị ban đầu bất kì s0 và p0 
 - Tính các giá trị bo,..,bn theo (1) 
 - Tính các giá trị co,...,cn theo (2) 
 - Tính ∆so và ∆po theo (3) và (4) 
 - Tính s1 = s0 + ∆so và p1 = po+ ∆po 
 - Lặp lại b−ớc 1 cho đến khi pi+1 = pi = p và si+1 = si = s 
 - Giải ph−ơng trình x2 - sx + p để tìm 2 nghiệm của đa thức 
- Bắt đầu quá trình trên cho đa thức Pn-2(x) 
Ví dụ : Tìm nghiệm của đa thức P4(x) = x
4 - 1.1x3 + 2.3x2 + 0.5x2 + 3.3. 
Với lần lặp ban đầu ta chọn s = -1 và p =1,nghĩa là tam thức có dạng x2 + x + 1 
 a0 a1 a2 a3 a4 
 1 -1.1 2.3 0.5 3.3 
 sbi -1 2.1 -3.4 0.8 
 -pbi-1 -1 2.1 -3.4 
 bi 1 -2.1 3.4 -0.8 = bn-1 0.7=bn 
 sbi -1.0 3.1 -5.5 
 -pbi-1 -1.0 3.1 
 ci 1 -3.1 5.5 -3.2 
 11.0
5.52.3
1.35.5
5.57.0
1.38.0
s =
−
−
−
−
=∆ 
06.0
5.52.3
1.35.5
7.02.3
8.05.5
p =
−
−
−=∆ 
s* = -1 + 0.11 = -0.89 
p* = 1 + 0.06 = 1.06 
Tiếp tục lặp lần 2 với s1 = s
* và p1 = p
* ta có : 
110
a0 a1 a2 a3 a4 
 1 -1.1 2.3 0.5 3.3 
 sbi -0.89 1.77 -2.68 0.06 
 -pbi-1 -1.06 2.11 -3.17 
 bi 1 -1.99 3.01 -0.07 = bn-1 0.17=bn 
 sbi -0.89 2.56 -4.01 
 -pbi-1 -1.0 3.1 
 ci 1 -2.88 4.51 -1.03 
01.0
51.403.1
88.251.4
5.57.0
88.207.0
s −=
−
−
−
−
=∆ 
04.0
51.403.1
88.251.4
17.003.1
07.051.4
p =
−
−
−−=∆ 
s* = -0.89 - 0.01 = -0.9 
p* = 1.06 + 0.04 = 1.1 
Nh− vậy P4(x) = (x2+0.9x+1.1)(x2 + 2x+3) 
Ch−ơng trình sau áp dụng lí thuyết vừa nêu để tìm nghiệm của đa thức. 
Ch−ơng trình 8-10 
//phuong phap Bairstow 
#include 
#include 
#include 
#include 
#define m 10 
void main() 
 { 
 float a[m],b[m],c[m]; 
 int i,n,v; 
 float s,e1,t,p,q,r,p1,q1; 
 clrscr(); 
 printf("Cho bac cua da thuc n = "); 
 scanf("%d",&n); 
 printf("Cho cac he so cua da thuc can tim nghiem\n"); 
 for (i=n;i>=0;i--) 
 { 
 printf("a[%d] = ",n-i); 
 scanf("%f",&a[i]); 
 } 
 printf("\n"); 
111
 e1=0.0001; 
 if (n<=2) 
 if (n==1) 
 { 
 printf("Nghiem cua he\n"); 
 printf("%.8f",(a[0]/(-a[1]))); 
 getch(); 
 exit(1); 
 } 
 do 
 { 
 v=0; 
 p=1; 
 q=-1; 
 b[n]=a[n]; 
 c[n]=a[n]; 
 do 
 { 
 b[n-1]=b[n]*p+a[n-1]; 
 c[n-1]=b[n-1]+b[n]*p; 
 for (i=n-2;i>=0;i--) 
 { 
 b[i]=b[i+2]*q+b[i+1]*p+a[i]; 
 c[i]=c[i+2]*q+c[i+1]*p+b[i]; 
 } 
 r=c[2]*c[2]-c[1]*c[3]; 
 p1=p-(b[1]*c[2]-b[0]*c[3])/r; 
 q1=q-(b[0]*c[2]-b[1]*c[1])/r; 
 if ((fabs(b[0])<e1)&&(fabs(b[1])<e1)) 
 goto tt; 
 v=v+1; 
 p=p1; 
 q=q1; 
 } 
 while (v<=40); 
 if(v>40) 
 { 
 printf("Khong hoi tu sau 40 lan lap"); 
 getch(); 
 exit(1); 
 } 
 tt:s=p1/2; 
 t=p1*p1+4*q1; 
 if(t<0) 
 { 
 printf("Nghiem phuc\n"); 
 printf("%.8f+%.8fj\n",s,(sqrt(-t)/2)); 
 printf("%.8f-%.8fj\n",s,(sqrt(-t)/2)); 
 printf("\n"); 
 } 
112
 else 
 { 
 printf("Nghiem thuc\n"); 
 printf("%.8f\n",(s+sqrt(t)/2)); 
 printf("%.8f\n",(s-sqrt(t)/2)); 
 printf("\n"); 
 } 
 for (i=2;i<=n;i++) 
 a[i-2]=b[i]; 
 n=n-2; 
 } 
 while ((n>2)&(r!=0.0)); 
 s=-a[1]/(2*a[2]); 
 t=a[1]*a[1]-4*a[2]*a[0]; 
 if (t<0) 
 { 
 printf("Nghiem phuc\n"); 
 printf("%.8f+%.8fj\n",s,(sqrt(-t)/(2*a[2]))); 
 printf("%.8f-%.8fj\n",s,(sqrt(-t)/(2*a[2]))); 
 printf("\n"); 
 } 
 else 
 { 
 printf("Nghiem thuc\n"); 
 printf("%.8f\n",(s-sqrt(t)/(2*a[2]))); 
 printf("%.8f\n",(s-sqrt(t)/(2*a[2]))); 
 printf("\n"); 
 } 
 getch(); 
 } 
 Dùng ch−ơng trình trên để xác định nghiệm của đa thức : 
 x6 - 2x5 - 4x4 + 13x3 - 24x2 + 18x - 4 = 0 
ta nhận đ−ợc các nghiệm : 
 x1 = 2.61903399 
 x2 = -2.73205081 
 x3 = 0.732050755 
 x4 = 0.381966055 
 x5 = 0.500011056 + i*1.3228881 
 x6 = 0.500011056 - i*1.3228881 
Đ11.Hệ ph−ơng trình phi tuyến 
 Ph−ơng pháp Newton có thể đ−ợc tổng quát hoá để giải hệ ph−ơng trình phi tuyến 
dạng : 
⎪⎪
⎪
⎩
⎪⎪
⎪
⎨
⎧
=
=
=
=
0)x,...,x,x,x(f
...............................
0)x,...,x,x,x(f
0)x,...,x,x,x(f
0)x,...,x,x,x(f
2321n
23213
23212
23211 
113
hay viết gọn hơn d−ới dạng : 
 F(X) = 0 
Trong đó : X = (x1,x2,x3,.....,xn) 
 Với một ph−ơng trình một biến,công thức Newton là : 
i i
i
i
x x
f x
f x+ = −1
( )
'( )
hay : f'(xi).∆x = -f(xi) 
với ∆x = xi+1 - xi 
Đối với hệ,công thức lặp là : 
 J(Xi)∆x = -F(Xi) 
Trong đó J(Xi) là toán tử Jacobi.Nó là một ma trận bậc n ( n - t−ơng ứng với số thành phần 
trong vectơ X) có dạng : 
i
n
n
n n n n
n
J x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
( ) .
.
..............
..............
..............
=
1
1
1
2
1
3
1
2
1
2
2
2
3
2
1 2 3
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
và ∆X = Xi+1 - Xi 
 Ph−ơng pháp Newton tuyến tính hoá hệ và nh− vậy với mỗi b−ớc lặp cần giải một hệ 
ph−ơng trình tuyến tính (mà biến là ∆xi) xác định bởi công thức lặp cho tới khi vectơ 
X(x1,x2,x3,.....,xn) gần với nghiệm. 
 D−ới đây là ch−ơng trình giải hệ ph−ơng trình phi tuyến 
 1
3
2
3
1 2 4
1 2 3 4
1
2
3
1 2 3 4
3 8 0
5 0
25 8 4 0
2 8 0
x x x x x
x x x x
x x
x x x x
− − =
+ + + − =
− + + =
− + =
−⎧
⎨
⎪⎪⎪
⎩
⎪⎪⎪
Ma trận đạo hàm riêng J(xi)là : 
1
2
2 4 2
2
1 4 1 2
1
2 3 2 3 2 3
3 3 3 3 0 3
1 1 1 1
25
0 8 0
2 2 2 1
1
2
x x x x x x x x
x
x
x x x x x x
− − −
−
−
−
−⎛
⎝
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎞
⎠
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
Ma trận này đ−ợc ch−ơng trình đọc vào nhờ thủ tục doc.Trong thủ tục này,các hệ số a[i,5] là 
các hàm fi(x).Vectơ nghiệm ban đầu đ−ợc chọn là { 0,-1,-1,1}T.Kết quả tính cho ta : x = 
{0.01328676,-1.94647929,-1.12499779,8.05819031 }T với độ chính xác 0.000001.Vectơ số 
d− r = { 0.00000536,-0.00000011,-0.00000001,-0.00000006}T. 
Ch−ơng trình 8-11 
//giai he pt phi tuyen 
#include 
#include 
114
#include 
#include 
#define n 4 
float a[n+1][n+2]; 
float x[n+1],y[n+1]; 
int i,j,k,l,z,r; 
float e,s,t; 
void main() 
 { 
 void doc(); 
 clrscr(); 
 printf("Cho cac gia tri nghiem ban dau\n"); 
 for (i=1;i<=n;i++) 
 { 
 printf("x[%d] = ",i); 
 scanf("%f",&x[i]); 
 } 
 e=1e-6; 
 z=30; 
 for (r=1;r<=z;r++) 
 { 
 doc(); 
 for (k=1;k<=n-1;k++) 
 { 
 s=0 ; 
 for (i=k;i<=n;i++) 
 { 
 t=fabs(a[i][k]); 
 if (s<=t) 
 { 
 s=t; 
 l=i; 
 } 
 } 
 for (j=k;j<=n+1;j++) 
 { 
 s=a[k][j]; 
 a[k][j]=a[l][j]; 
 a[l][j]=s; 
 } 
 if (a[1][1]==0) 
 { 
 printf("Cac phan tu duong cheo cua ma tran bang khong"); 
 getch(); 
 exit(1); 
 } 
 else 
 { 
115
 if (fabs(a[k][k]/a[1][1])<(1e-08)) 
 { 
 printf("Ma tran suy bien"); 
 goto mot; 
 } 
 } 
 for (i=k+1;i<=n;i++) 
 { 
 if (a[k][k]==0) 
 { 
 printf("Cac phan tu duong cheo cua ma tran bang 
khong\n"); 
 goto mot; 
 } 
 s=a[i][k]/a[k][k]; 
 a[i][k]=0; 
 for (j=k+1;j<=n+1;j++) 
 a[i][j]=a[i][j]-s*a[k][j]; 
 } 
 y[n]=a[n][n+1]/a[n][n]; 
 for (i=n-1;i>=1;i--) 
 { 
 s=a[i][n+1]; 
 for (j=i+1;j<=n;j++) 
 s=s-a[i][j]*y[j]; 
 if (a[i][i]==0) 
 { 
 printf("Cac phan tu duong cheo cua ma tran bang 
khong\n"); 
 goto mot; 
 } 
 y[i]=s/a[i][i]; 
 } 
 } 
 if (r!=1) 
 for (i=1;i<=n;i++) 
 { 
 if (fabs(y[i])<e*fabs(x[i])) 
 goto ba; 
 } 
 for (i=1;i<=n;i++) 
 x[i]=x[i]-y[i]; 
 printf("\n"); 
 } 
 printf("Khong hoi tu sau %d lan lap\n",z); 
 goto mot; 
 clrscr(); 
 ba:printf("Vec to nghiem\n"); 
 for (i=1;i<=n;i++) 
 printf("%.5f\n",(x[i]-y[i])); 
116
 printf("\n"); 
 printf("Do chinh xac cua nghiem la %.5f: \n", e); 
 printf("\n"); 
 printf("Vec to tri so du :\n"); 
 for (i=1;i<=n;i++) 
 printf("%.5f\n",(a[i][n+1])); 
 mot:printf("\n"); 
 getch(); 
 } 
void doc() 
 { 
 a[1][1]=3*x[1]*x[1]-3*x[2]*x[4]; 
 a[1][2]=-3*x[2]*x[2]-3*x[1]*x[4]; 
 a[1][3]=0; 
 a[1][4]=-3*x[1]*x[2]; 
 a[1][5]=x[1]*x[1]*x[1]-x[2]*x[2]*x[2]-3*x[1]*x[2]*x[4]-8; 
 a[2][1]=1; 
 a[2][2]=1; 
 a[2][3]=1; 
 a[2][4]=1; 
 a[2][5]=x[1]+x[2]+x[3]+x[4]-5; 
 a[3][1]=-x[1]/sqrt(25-x[1]*x[1]); 
 a[3][2]=0; 
 a[3][3]=8; 
 a[3][4]=0; 
 a[3][5]=sqrt(25-x[1]*x[1])+8*x[3]+4; 
 a[4][1]=2*x[2]*x[3]; 
 a[4][2]=2*x[1]*x[3]; 
 a[4][3]=2*x[1]*x[2]; 
 a[4][4]=-1; 
 a[4][5]=2*x[1]*x[2]*x[3]-x[4]+8; 
 }