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

87 Ch−ơng 8 : Giải gần đúng ph−ơng trình đại số và siêu việt Đ1.Khái niệm chung Nếu ph−ơng trình đại số hay siêu việt khá phức tạp thì ít khi tìm đ−ợc nghiệm đúng.Bởi vậy việc tìm nghiệm gần đúng và −ớc l−ợng sai số là rất cần thiết. Ta xét ph−ơng trình : f(x) = 0 (1) với f(x) là hàm cho tr−ớc của biến x.Chúng ta cần tìm giá trị gần đúng của nghiệm của ph−ơng trình này. Quá trình giải th−ờng chia làm hai b−ớc: b−ớc sơ bộ và b−ớc kiện toàn nghiệm. B−ớc giải sơ bộ có 3 nhiệm vụ:vây nghiệm, tách nghiệm và thu hẹp khoảng chứa nghiệm. Vây nghiệm là tìm xem các nghiệm của ph−ơng trình có thể nằm trên những đoạn nào của trục x.Tách nghiệm là tìm các khoảng chứa nghiệm soa cho trong mỗi khoảng chỉ có đúng một nghiệm.Thu hẹp khoảng chứa nghiệm là làm cho khoảng chứa nghiệm càng nhỏ càng tốt.Sau b−ớc sơ bộ ta có khoảng chứa nghiệm đủ nhỏ. B−ớc kiện toàn nghiệm tìm các nghiệm gần đúng theo yêu cầu đặt ra. Có rất nhiều ph−ơng pháp xác định nghiệm của (1).Sau đây chúng ta xét từng ph−ơng pháp. Đ2.Ph−ơng pháp lặp đơn Giả sử ph−ơng trình (1) đ−ợc đ−a về dạng t−ơng đ−ơng : x = g(x) (2) từ giá trị xo nào đó gọi là giá trị lặp đầu tiên ta lập dãy xấp xỉ bằng công thức : xn = g(xn-1) (3) với n = 1,2,.... Hàm g(x) đ−ợc gọi là hàm lặp.Nếu dãy xn → α khi n →∝ thì ta nói phép lặp (3) hội tụ. x1 xo xo x1 Hình a Hình b Ta có định lí:Xét ph−ơng pháp lặp (3),giả sử : - [a,b] là khoảng phân li nghiệm α của ph−ơng trình (1) tức là của (2) - mọi xn tính theo (3) đều thuộc [a,b] - g(x) có đạo hàm thoả mãn : 88 bxa,1q)x(g <<<≤′ (4) trong đó q là một hằng số thì ph−ơng pháp lặp (3) hội tụ Ta có thể minh hoạ phép lặp trên bằng hình vẽ a và b. Cách đ−a ph−ơng trình f(x) = 0 về dạng x = g(x) đ−ợc thực hiện nh− sau:ta thấy f(x) = 0 có thể biến đổi thành x = x + λf(x) với λ ≠ 0.Sau đó đặt x + λf(x) = g(x) sao cho điều kiện (4) đ−ợc thoả mãn. Ví dụ:xét ph−ơng trình x3 + x - 1000 = 0 Sau b−ớc giải sơ bộ ta có nghiệm x1 ∈ ( 9,10 ) Nếu đ−a ph−ơng trình về dạng: x = 1000 - x3 = g(x) thì dễ thấy | g'(x) | > 1 trong khoảng ( 9,10 ) nên không thoả mãn điều kiện (4) Chúng ta đ−a ph−ơng trình về dạng 3 x1000x −= thì ta thấy điều kiện (4) đ−ợc thoả mãn.Xây dựng dãy xấp xỉ 3 n1n x1000x −=+ với xo chọn bất kì trong ( 9,10 ) Trên cơ sở ph−ơng pháp này chúng ta có các ch−ơng trình tính toán sau: Ch−ơng trình giải ph−ơng trình exp((1/3)*ln(1000-x)) với số lần lặp cho tr−ớc Ch−ơng trình 8-1 //lap don #include #include #include void main() { int i,n; float x,x0; float f(float); clrscr(); printf("Cho so lan lap n = "); scanf("%d",&n); printf("Cho gia tri ban dau cua nghiem x0 = "); scanf("%f",&x0); x=x0; for (i=1;i<=n;i++) x=f(x); printf("Nghiem cua phuong trinh la :%.4f",x); getch(); } float f(float x) { float a=exp((1./3.)*log(1000-x)); return(a); 89 } và ch−ơng trình giải bài toán bằng ph−ơng pháp lặp với sai số cho tr−ớc Ch−ơng trình 8-2 //lap don #include #include #include void main() { int i; float epsi,x,x0,y; float f(float); clrscr(); printf("Cho sai so epsilon = "); scanf("%f",&epsi); printf("Cho gia tri ban dau cua nghiem x0 = "); scanf("%f",&x0); x=x0; y=f(x); if (abs(y-x)>epsi) { x=y; y=f(x); } printf("Nghiem cua phuong trinh la %.6f",y); getch(); } float f(float x) { float a=exp((1./3.)*log(1000-x)); return(a); } Cho giá trị đầu xo = 1.Kết quả tính toán x = 9.966555 Đ3.Ph−ơng pháp chia đôi cung 90 Giả sử cho ph−ơng trình f(x) = 0 với f(x) liên tục trên đoạn [a,b] và f(a).f(b) < 0.Chia đoạn [a,b] thành 2 phần bởi chính điểm chia (a + b)/2. 1.Nếu f((a+b)/2) = 0 thì ξ = (a+b)/2 2.Nếu f((a+b)/2) ≠ 0 thì chọn [ a,(a + b)/2 ] hay [ (a + b)/2,b ] mà giá trị hàm hai đầu trái dấu và kí hiệu là [a1,b1].Đối với [a1,b1] ta lại tiến hành nh− [a,b] Cách làm trên đ−ợc mô tả trong ch−ơng trình sau dùng để tìm nghiệm của ph−ơng trình : x4 + 2x3 - x - 1 = 0 trên đoạn [0,1] Ch−ơng trình 8-3 //chia doi cung #include #include #include #define epsi 0.00001 void main() { float x0,x1,x2; float y0,y1,y2; float f(float); int maxlap,demlap; clrscr(); printf("Tim nghiem cua phuong trinh phi tuyen"); printf("\nbang cach chia doi cung\n"); printf("Cho cac gia tri x0,x1,maxlap\n"); printf("Cho gia tri x0 = "); scanf("%f",&x0); printf("Cho gia tri x1 = "); scanf("%f",&x1); printf("Cho so lan lap maxlap = "); scanf("%d",&maxlap); y0=f(x0); y1=f(x1); if ((y0*y1)>0) { printf("Nghiem khong nam trong doan x0 - x1\n"); printf(" x0 = %.2f\n",x0); printf(" x1 = %.2f\n",x1); printf(" f(x0) = %.2f\n",y0); printf(" f(x1) = %.2f\n",y1); } demlap=0; do { y x a b ξ b1 91 x2=(x0+x1)/2; y2=f(x2); y0=f(x0); if (y0*y2>0) x0=x2; else x1=x2; demlap=demlap+1; } while(((abs((y2-y0))>epsi)||(demlap<maxlap))); if (demlap>maxlap) { printf("Phep lap khong hoi tu sau %d lan lap ",maxlap); printf(" x0 = %.2f\n",x0); printf(" x1 = %.2f\n",x1); printf(" f(x2) = %.2f\n",y2); } else { printf("Phep lap hoi tu sau %d lan lap\n",demlap); printf("Nghiem x = %.2f",x2); } getch(); } float f(float x) { float a=x*x*x*x+2*x*x*x-x-1 ; return(a); } Kết quả tính cho nghiệm:x = 0.87 Đ4.Ph−ơng pháp dây cung Giả sử f(x) liên tục trên trên đoạn [a,b] và f(a).f(b) < 0.Cần tìm nghiệm của f(x) = 0.Để xác định ta xem f(a) 0.Khi đó thay vì chia đôi đoạn [a,b] ta chia [a,b] theo tỉ lệ -f(a)/f(b).Điều đó cho ta nghiệm gần đúng : x1 = a + h1 Trong đó 1h f a f a f b b a= − − + − ( ) ( ) ( ) ( ) Tiếp theo dùng cách đó với đoạn [ a,x1] hay [ x1,b] mà hai đầu hàm nhận giá trị trái dấu ta đ−ợc nghiệm gần đúng x2 v.v. Về mặ hình học,ph−ơng pháp này có nghĩa là kẻ dây cung của đ−ờng cong f(x) qua hai điểm A[a,f(a)] và B[b,f(b)].Thật vậy ph−ơng trình dây cung AB có dạng : 92 )a(f)b(f )a(fy ab ax − −=− − Cho x = x1 y = 0 ta có )ab()a(f)b(f )a(f ax1 −−−= Trên cơ sở của ph−ơng pháp ta có ch−ơng trình tính nghiệm của ph−ơng trình x4 + 2x3 - x - 1 = 0 trên đoạn [0,1] Ch−ơng trình 8-4 //phuong phap day cung #include #include #include #define epsi 0.00001 void main() { float a,b,fa,fb,dx,x; float f(float); clrscr(); printf("Tim nghiem cua phuong trinh phi tuyen\n"); printf("bang phuong phap day cung\n"); printf("Cho cac gia tri a,b\n"); printf("Cho gia tri cua a = "); scanf("%f",&a); printf("Cho gia tri cua b = "); scanf("%f",&b); fa=f(a); fb=f(b); dx=fa*(b-a)/(fa-fb); while (fabs(dx)>epsi) { x=a+dx; fa=f(x); if((fa*fb)<=0) a=x; else b=x; fa=f(a); fb=f(b); dx=fa*(b-a)/(fa-fb); } B b a x1 ξ A 93 printf("Nghiem x = %.3f",x); getch(); } float f(float x) { float e=x*x*x*x+2*x*x*x-x-1; return(e); } Kết quả tính cho nghiệm:x = 0.876 Đ5.Ph−ơng pháp lặp Newton Ph−ơng pháp lặp Newton (còn gọi là ph−ơng pháp tiếp tuyến) đ−ợc dùng nhiều vì nó hội tụ nhanh.Giả sử f(x) có nghiệm là ξ đã đ−ợc tách trên đoạn [a,b] đồng thời f'(x) và f"(x) liên tục và giữ nguyên dấu trên đoạn [a,b].Khi đã tìm đ−ợc xấp xỉ nào đó xn ∈ [a,b] ta có thể kiện toàn nó theo ph−ơng pháp Newton.Từ mút B ta vẽ tiếp tuyến với đ−ờng cong.Ph−ơng trình đ−ờng tiếp tuyến là )xx)(b(f)x(fy 00 −′=− Tiếp tuyến này cắt trục hoành tại điểm có y=0,nghĩa là : )xx)(b(f)x(f 010 −′=− hay : )x(f )x(f xx 0 0 01 ′−= Từ x1 ta lại tiếp tục vẽ tiếp tuyến với đ−ờng cong thì giao điểm xi sẽ tiến tới nghiệm của ph−ơng trình. Việc chọn điểm ban đầu xo rất quan trọng.Trên hình vẽ trên ta thấy nếu chọn điểm ban đầu xo = a thì tiếp tuyến sẽ cắt trục tại một điểm nằm ngoài đoạn [a,b].Chọn xo = b sẽ thuận lợi cho việc tính toán.Chúng ta có định lí : Nếu f(a).f(b) < 0 ; f(x) và f"(x) khác không và giữ nguyên dấu xác định khi x ∈ [a,b] thì xuất phát từ xo∈ [a,b] thoả mãn điều kiện f(xo).f″(xo) > 0 có thể tính theo ph−ơng pháp Newton nghiệm ξ duy nhất với độ chính xác tuỳ ý. Khi dùng ph−ơng pháp Newton cần lấy xo là đầu mút của đoạn [a,b] để tại đó f(xo).f"(xo) > 0.áp dụng lí thuyết trên chúng ta xây dựng ch−ơng trình tính sau: Ch−ơng trình 8-5 //phuong phap Newton #include #include a b=xox1 94 #include #include #define n 50 #define epsi 1e-5 void main() { float t,x0; float x[n]; int i; float f(float); float daoham(float); clrscr(); printf("Tim nghiem cua phuong trinh phi tuyen\n"); printf("bang phuong phap lap Newton\n"); printf("Cho gia tri cua x0 = "); scanf("%f",&x0); i=1; x[i]=x0; do { x[i+1] = x[i]-f(x[i])/daoham(x[i]); t = fabs(x[i+1]-x[i]); x[i]=x[i+1]; i=i+1; if (i>100) { printf("Bai toan khong hoi tu\n"); getch(); exit(1); } else ; } while (t>=epsi); printf("Nghiem x = %.5f",x[i]); getch(); } float f(float x) { float a=x*x-x-2; return(a); } float daoham(float x) { float d=2*x-1; return(d); 95 } Ch−ơng trình này đ−ợc dùng xác định nghiệm của hàm đã đ−ợc định nghĩa trong function.Trong tr−ờng hợp này ph−ơng trình đó là:x2 - x -1 = 0.Kết quả tính với giá trị đầu xo = 0 cho nghiệm x = 2. Đ6.Ph−ơng pháp Muller Trong ph−ơng pháp dây cung khi tìm nghiệm trong đoạn [a,b] ta xấp xỉ hàm bằng một đ−ờng thẳng.Tuy nhiên để giảm l−ợng tính toán và để nghiệm hội tụ nhanh hơn ta có thể dùng ph−ơng pháp Muller.Nội dung của ph−ơng pháp này là thay hàm trong đoạn [a,b] bằng một đ−ờng cong bậc 2 mà ta hoàn toàn có thể tìm nghiêm chín xác của nó.Gọi các điểm đó có hoành độ lần l−ợt là a = x2,b = x1 và ta chọn thêm một điểm x0 nằm trong đoạn [x2,x1].Gọi h1 = x1 - x0 h2 = x0 - x2 v = x - x0 f(x0) = f0 f(x1) = f1 f(x2) = f2 1 2 h h=γ Qua 3 điểm này ta có một đ−ờng parabol : y = av2 + bv + c Ta tìm các hệ số a,b,c từ các giá trị đã biết v: 22 2 222 11 2 111 0 2 0 fcbhah)xx(hv fcbhah)xx(hv fc)0(b)0(a)xx(0v =++== =++== =++== Từ đó ta có : 0 1 2 101 2 1 201 fc h ahff b )1(h f)1(ff a = −−= + ++−= γγ γγ Sau đó ta tìm nghiệm của ph−ơng trình av2 + bv + c = 0 và có : ac4bb c2 xn 202,1 −± −= Tiếp đó ta chọn nghiệm gần x0 nhất làm một trong 3 điểm để tính xấp xỉ mới.Các điểm này đ−ợc chọn gần nhau nhất. Tiếp tục quá trình tính đến khi đạt độ chính xác yêu cầu thì dừng lại. Ví dụ:Tìm nghiệm của hàm f(x) = sin(x) - x/2 trong đoạn [1.8,2.2].Ta chọn x0 = 2 Ta có : x0 = 2 f(x0) = -0.0907 h1 = 0.2 x1 = 2.2 f(x1) = -0.2915 h2 = 0.2 x2 = 1.8 f(x2) = 0.07385 γ = 1 Vậy thì : x1,f1 x0,f0 x2,f2 av 2+bv+c f(x) h1 h2 96 0907.0c 91338.0 2.0 2.0)45312.0()097.0(2915.0 b 45312.0 )11(2.01 07385.0)11()0907.0()2915.0(1 a 2 2 −= −=ì−−−−−= −=+ìì ++ì−−−ì= Ta có nghiệm gần x0 nhất là : 89526.1 )0907.0()45312.0(4)91338.0(91338.0 )0907.0(2 0.2n 21 = −ì−ì−−−− −ì−= Với lần lặp thứ hai ta có : x0 = 1.89526 f(x0) = 1.9184ì10-4 h1 = 0.10474 x1 = 2.0 f(x1) = -0.0907 h2 = 0.09526 x2 = 1.8 f(x2) = 0.07385 γ = 0.9095 Vậy thì : 4 24 2 4 109184.1c 81826.0 10474.0 10474.0)4728.0(109184.10907.0 b 4728.0 9095.110474.09095.0 07385.09095.1)109184.1()0907.0(9095.0 a − − − ì= −=ì−−ì−−= −=ìì +ìì−−ì= Ta có nghiệm gần x0 nhất là : 89594.1 109184.1)4728.0(4)81826.0(81826.0 109184.12 89526.1n 42 4 1 =ìì−ì−−− ìì−= − − Ta có thể lấy n1 = 1.895494 làm nghiệm của bài toán Ch−ơng trình giải bài toán bằng ph−ơng pháp Muller nh− sau: Ch−ơng trình 8-6 //phuong phap Muller #include #include #include #include void main() { float x0,x1,x2,h1,h2,eps; float a,b,c,gamma,n1,n2,xr; int dem; float f(float); clrscr(); printf("PHUONG PHAP MULLER\n"); printf("\n"); printf("Cho khoang can tim nghiem [a,b]\n"); printf("Cho gia tri duoi a = "); scanf("%f",&x2); 97 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; 98 } } 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 : aox n + a1x n-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 α 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 α ta chia đa thức Pn(x) cho (x - α) 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 ph−ơng trình sai phân ϕ có dạng nh− sau : ϕ = 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) = a0x n + a1x n-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à xy k ii = Nếu yi là các nghiệm của ph−ơng trình sai phân là tuyến tính (1),thì xcxcxcy knn k 22 k 11k ....+++= (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 là sao cho : | x1| ≥ | x2 | ≥...| xn| Vậy thì k k ky c x c c x x = + +1 1 1 2 2 1 1[ ( ) ]... và ]) x x( c c1[xcy ... 1k 1 2 2 11k 111k ++= +++ 99 do đó : ]) x x( c c1[ ]) x x( c c1[ xy y ... ... k 1 2 1 2 1k 1 2 1 2 1 k 1k + + = + + + + do x1 > x2 nên: ∞→→ kkhi0) x x() x x( ......., k 1 2k 1 2 vậy thì : ∞→∞→+ kkhi y y k 1k Nghĩa là : y y limx k 1k k 1 + ∞→ = 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 : )yaya( a 1y kn1nk1 o nk ... ++−= −++ (4) Ví dụ:Tính nghiệm của đa thức Pn(x) = P3(x) = x 3 - 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 8-7 //phuong phap Bernoulli #include #include 100 #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; 101 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 )x(P )x(P xx n i in i1i ′−=+ (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ụ : a a x 1n n o − −= và tính tiếp các giá trị sau : )x(P )x(P xx n o on o1 ′−= )x(P )x(P xx n 1 1n 12 ′−= 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) = box n-1 + b1x n-2+p3x n-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)( box n-1 + b1x n-2+b3x n-3 +..+ bn-2x + bn-1 ) + bn = aox n + a1x n-1 + a2x n-2 +...+ an-1x + an (5) Từ (5) rút ra : bo = ao b1 = a1 + boxi (6) b2 = a2 + b1xi ...... bk = ak + bk-1xi ....... 102 bn = an + bn-1xi = Pn(xi) Đạo hàm (3) ta đ−ợc : )x(P)x(P)xx()x(P 1n1nin −− +′−=′ và )x(P)x(P i1nin −=′ (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 | < ε hay Pn(xi+1) ≈ 0 nên α1≈ xi+1 là một nghiệm của đa thức. Phép chia Pn(x) cho (x - α1) 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 = α1.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) = x 3 - 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) = b0x 2 + b1x + b2 = 13.75 6.3 75.13 375.1 5.3 )x(P )x(P xx 0n 0n 01 =+=′−= 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) = b0x 2 + b1x + b2 = 15.68 606.3 68.15 096.0 6.3 )x(P )x(P xx 1n 1n 12 =+=′−= 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 8-8 //phuong phap Birge-Viette #include #include #include #define max 20 103 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; } 104 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) 105 Do (6) đúng với mọi n nên cho n = 1 , 2 , 3 , . . . ta có : | x2 - x1 | ≤ q | x1 - xo | | x3 - x2 | ≤ q | x2 - x1 | . . . . . . . . . . . . . . . . . . . | xn+1 - xn | ≤ q | xn - xn-1 | Điều này có nghĩa là dãy xi+1 - xi , một cách gần đúng,là một cấp số nhân . Ta cũng coi rằng dãy xn - y với y là nghiệm đúng của (1) , gần đúng nh− một cấp số nhân có công sai q . Nh− vậy : 1q yx yx n 1n <=− −+ (7) hay : )yx(qyx n1n −=−+ (8) T−ơng tự ta có : )yx(qyx 1n2n −=− ++ (9) Từ (8) và (9) ta có : n1n 1n2n xx xx q − −= + ++ (10) Thay giá trị của q vừa tính ở (10) vào biểu thức của q ở trên ta có : ( ) 1n1nn 2 1nn n xx2x xx xy ++ + +− −−= (11) Công thức (11) đ−ợc gọi là công thức ngoại suy Adam.Nh− vậy theo (11) tr−ớc hết ta dùng ph−ơng pháp lặp để tính giá trị gần đúng xn+2,xn+1,xn của nghiệm và sau đó theo (11) tìm đ−ợc nghiệm với sai số nhỏ hơn. Để làm ví dụ chúng ta xét ph−ơng trình : lnx - x2 + 3 = 0 Ta đ−a về dạng lặp : 3)xln(x += 3xlnx2 1 )x(f +=′ Phép lặp hội tụ trong đoạn [0.3,∝].Ta cho x1 = 1 thì tính đ−ợc : x2 = 1,7320508076 x3 = 1.883960229 x4 = 1.90614167 y = 1.909934347 Để giảm sai số ta có thể lặp nhiều lần Ch−ơng trình 8-9 //phuong phap Aitken #include #include #include #define m 5 void main() { float x[m]; 106 float epsi,n,y; int i,z; float f(float); clrscr(); printf("Cho tri so ban dau x[1] = "); scanf("%f",&x[1]); printf("Cho tri so sai so epsilon = "); scanf("%f",&epsi); printf("\n"); printf( "Ngoai suy Aitken cua ham\n"); z=0; while (z<=20) { for (i=2;i<=4;i++) x[i]=f(x[i-1]); n=x[4]-2*x[3]+x[2]; if ((fabs(n)<1e-09)||(fabs(x[1]-x[2])<epsi*fabs(x[1]))) z=20; else { y=x[2]-(x[3]-x[2])*(x[3]-x[2])/n; if (z>20) printf("Khong hoi tu sau hai muoi lan lap\n"); x[1]=y; } z=z+1; } printf("Nghiem cua phuong trinh y = %.6f",y); getch(); } float f(float x) { float s=sqrt(log(x)+3); return(s); } Với giá trị ban đầu là 1 và sai số là 1e-8,ch−ơng trình cho kết quả y = 1.9096975944 Đ10.Ph−ơng pháp Bairstow Nguyên tắc của ph−ơng pháp Bairstow là trích từ đa thức Pn(x) một tam thức Q2(x) = x2 - sx + p mà ta có thể tính nghiệm thực hay nghiệm phức của nó một cách đơn giản bằng các ph−ơng pháp đã biết. Việc chia đa thức Pn(x) cho tam thức Q2(x) đ−a tới kết quả : Pn(x) = Q2(x).Pn-2(x) + R1(x) với Pn(x) = aox n + a1x n-1 + a2x n-2 +...+ an Q2(x) = x 2 - sx + p 107 Pn-2(x) = box n-2 + b1x n-3 + b2x n-4 +...+ bn-2 R1(x) = αx + β Để có đ−ợc một th−ơng đúng,cần tìm các giá trị của s và p sao cho R1(x) = 0 (nghĩa là α và β triệt tiêu).Với s và p đã cho,các hệ số b của đa thức Pn-2(x) và các hệ số α và β đ−ợc tính bằng ph−ơng pháp truy hồi.Các công thức nhận đ−ợc khi khai triển biểu thức Pn(x) = Q2(x).Pn-2(x) + R1(x) và sắp xếp lại các số hạng cùng bậc : aox n + a1x n-1 + a2x n-2 +...+ an = (x 2 - sx + p)( box n-2 + b1x n-3 + b2x n-4 +...+ bn-2) Số hạng bậc Hệ số của Pn(x) Hệ số của Q2(x).Pn-2(x) xn ao bo xn-1 a1 b1 - sbo xn-2 a2 b2 - sb1 + pbo ...... ...... ..... xn-k ak bk - sbk-1 + pbk-2 x an-1 α - sbn-2 + pbn-3 xo an β + pbn-2 Nh− vậy : bo = a o (1) b1 = a1 + sbo b2 = a2 + sb1 - pbo .................. bk = ak + sbk-1 - pbk-2 α = an-1 + sbn-2 - pbn-3 β = an - pbn-2 Chúng ta nhận thấy rằng α đ−ợ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 = α Hệ số bn là : bn = an + sbn-1 - pbn-2 = sbn-1 + β và cuối cùng : R1(x) = αx + β = bn-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 x 2 - 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: ⎩⎨ ⎧ ==0)p,s(g 0)p,s(f 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) 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; }