Giáo trình Matlab cơ bản - Chương 5: Các phương trình phi tuyến

i) 
với   Xi = { ri, si}T và   Xi+1 = { ri+1, si+1}T 
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
i i
i
i i
f(r ,s )
F(X )
g(r ,s )
  i
f f
r sJ(X )
g g
r s
∂ ∂⎛ ⎞⎜ ⎟∂ ∂= ⎜ ⎟∂ ∂⎜ ⎟⎜ ⎟∂ ∂⎝ ⎠
Quan hệ  :  J(Xi)∆X  =  ‐F(Xi) với ∆X  =  {ri+1  ‐  ri,si+1  ‐  si}T  tương  ứng với một hệ 
phương trình tuyến tính hai ẩn số ∆r = ri+1 ‐ ri  và  ∆s = si+1 ‐ si : 
i i
i i
f fr s f(r ,s )
r s
g gr s g(r ,s )
r s
∂ ∂⎧ ∆ + ∆ = −⎪⎪∂ ∂⎨∂ ∂⎪ ∆ + ∆ = −⎪ ∂ ∂⎩
Theo công thức Cramer ta có : 
302
g ff g
s sr
∂ ∂− +∂ ∂∆ = δ  
gfg f
r rs
∂∂− +∂ ∂∆ = δ  
  g gf f
r s s r
∂ ∂∂ ∂δ= −∂ ∂ ∂ ∂  
Để dùng  được  công  thức này  ta cần  tính  được các  đạo hàm  f
r
∂
∂ ,
f
s
∂
∂ ,
g
r
∂
∂ ,
g
s
∂
∂ . 
Các đạo hàm này được tính theo công thức truy hồi. 
Do b1 = a1 nên  
  1 1b a 0
r r
∂ ∂= =∂ ∂      
1 1b a 0
s s
∂ ∂= =∂ ∂       
b2 = a2 + rb1  nên   
2 1
1 1
b br b b
r r
∂ ∂= + =∂ ∂    
2 2 1b a br 0
s s s
∂ ∂ ∂= + =∂ ∂ ∂    
b3 = a3 + rb2 + sb1 nên      
3 3 2 1b a (rb ) (sb )
r r r r
∂ ∂ ∂ ∂= + +∂ ∂ ∂ ∂  
Mặt khác :   
3a 0
r
∂ =∂  
2 2
2
(rb ) (b )r b
r r
∂ ∂= +∂ ∂     
1(rb ) 0
r
∂ =∂        
nên:   3 2 1
b b rb
r
∂ = +∂      
3 2 1
1 1
b b br s b b
s s s
∂ ∂ ∂= + + =∂ ∂ ∂      
b4 = a4 + rb3 + sb2 nên:   
4 3 2
3
b b bb r s
r r r
∂ ∂ ∂= + +∂ ∂ ∂    
4 3 2
2
b b br s b
s s s
∂ ∂ ∂= + +∂ ∂ ∂  
. . . . .  
n n 1 n 2
n 1
b b br s b
r r r
− −
−
∂ ∂ ∂= + +∂ ∂ ∂  
n n 1 n 2
n 2
b b br s b
s r r
− −
−
∂ ∂ ∂= + +∂ ∂ ∂  
Nếu chúng ta đặt: 
  k k 1
b c
r −
∂ =∂        
thì :   
c1 = b1                    (2) 
  c2 = b2 + rc1 
  c3 = b3 + rc2 + sc1 
303
  . . . . . . . . . .  
  cn = bn + rcn‐1 + scn‐2 
Như vậy ta có: 
n n 1 n 1 n 2
2
n 1 n n 2
b c b cr
c c c
− + −
− −
− +∆ = −                 (3) 
n 1 n 1 n n
2
n 1 n n 2
b c b cs
c c c
+ −
− −
− +∆ = −                 (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ị b1,..,bn+1 theo (1) 
  ‐ Tính các giá trị c1,...,cn theo (2) 
  ‐ Tính ∆ro và ∆so theo (3) và (4) 
  ‐ Tính s1 = r0 + ∆ro và s1 = so+ ∆so 
  ‐ Lặp lại bước 1 cho đến khi ri+1 = ri = r và si+1 = si = s 
  ‐ Giải phương trình x2 ‐ rx ‐ s để 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) 
Ta xây dựng hàm bairstow() để thực hiện thuật toán tìm r, s 
function [r,s] = bairstow(p, r0, s0, maxiter) 
% tim da thuc bac 2 dang x^2 ‐ rx ‐ s 
%vao ‐p la da thuc can tim nghiem 
%    ‐r0, s0 gia tri ban dau 
%    ‐maxiter so lan lap max 
%ra  ‐r, s 
%cu phap [r,s] = bairstow(p, r0, s0 ,maxiter) 
n = length(p) ‐ 1; 
c = zeros(n); 
b = zeros(n); 
j = 0; 
while j < maxiter 
    b(1) = p(1); 
    b(2) = p(2) + r0*b(1); 
    for k = 3:(n+1) 
304
        b(k) = p(k) + r0*b(k‐1) + s0*b(k‐2); 
    end 
    c(1) = b(1); 
    c(2) = b(2) + r0*c(1); 
    for k = 3:(n) 
        c(k) = b(k) + r0*c(k‐1) + s0*c(k‐2); 
    end 
    d0 = det([c(n‐1), c(n‐2); c(n), c(n‐1)]); 
    d1 = det([‐b(n), c(n‐2); ‐b(n+1), c(n‐1)]); 
    d2 = det([c(n‐1), ‐b(n); c(n) ‐b(n+1)]); 
    r = r0 + d1/d0; 
    s = s0 + d2/d0; 
    if ((abs(d1/d0))&(abs(d2/d0)))<eps 
        break; 
    end 
    r0 = r; 
    s0 = s; 
    j = j + 1; 
end 
Để tìm nghiệm của đa thức P4(x) = x4 ‐ 1.1x3 + 2.3x2 ‐ 0.5x + 3.3 ta dùng chương 
trình ctbairstow.m 
 clear all, clc 
 p = [1  ‐1.1  2.3  0.5  3.3]; 
 m = length(p); 
 s0 = ‐1,; 
 r0 = ‐1; 
 fprintf(ʹNghiem cua da thuc:\nʹ); 
 while m > 3 
    [r, s] = bairstow(p,r0,s0,50); 
    m = m ‐ 2; 
    x1 = (r + sqrt(r^2+4*s))/2; 
    x2 = (r ‐ sqrt(r^2+4*s))/2; 
    fprintf(ʹ%s\nʹ,num2str(x1)); 
    fprintf(ʹ%s\nʹ,num2str(x2)); 
305
    p = deconv(p,[1 ‐r ‐s]); 
 end 
 if length(p) == 3 
    x1 = (‐p(2) + sqrt(p(2)^2‐4*p(3)))/2; 
    x2 = (‐p(2) ‐ sqrt(p(2)^2‐4*p(3)))/2; 
    fprintf(ʹ%s\nʹ,num2str(x1));; 
    fprintf(ʹ%s\nʹ,num2str(x2));;  
 else 
    x1 = ‐p(2)/p(1); 
    fprintf(ʹ%f\nʹ,x1);; 
 end 
§19. PHƯƠNG PHÁP LOBACHEVSKY ‐ GRAEFFE 
  Phương  pháp  này  đã  được  Karl  Heinrich  Gräffe,  Germinal  Pierre 
Dandelin và Nikolai Ivanovich Lobachevsky đưa ra. Nó có một nhược điểm 
là các kết quả trung gian có trị số rất lớn. 
Xét phương trình : 
     P(x) = a0xn + a1xn‐1 + ⋅ ⋅ ⋅ + an = 0          (1) 
chỉ có nghiệm thực với giá trị tuyệt đối khác nhau. Ta đánh số các nghiệm đó 
theo thứ tự giá trị tuyệt đối giảm : 
    |x1| > |x2| > ⋅ ⋅ ⋅ > |xn|              (2) 
Dựa vào (1) ta xây dựng một phương trình mới : 
    Q(x) = c0xn + c1xn‐1 + ⋅ ⋅ ⋅ + cn = 0  ( c0 ≠ 0)      (3) 
có nghiệm là   m m m1 2 nx , x ,..., x− − − . Sau đó ta viết lại phương trình trên: 
            m m m0 1 2 nQ(x) c (x x )(x x ) (x x )= + + ⋅ ⋅ ⋅ +          (4) 
So sánh (3) và (4) ta có: 
m m m 1
1 2 n
0
m m m m m m 2
1 2 2 3 n 1 n
0
m m m m n
1 2 n 1 n
0
cx x x
c
cx x x x x x
c
cx x x x
c
−
−
⎧ + + ⋅ ⋅ ⋅ + =⎪⎪⎪ + + ⋅ ⋅ ⋅ + =⎪⎨⎪⋅ ⋅ ⋅⎪⎪ ⋅ ⋅ ⋅ =⎪⎩
          (5) 
Vì có giả thiết (2) nên khi m tăng thì ở vế trái của các đẳng thức (5) số hạng 
đầu trội lên, lấn át các số hạng sau và với m đủ lớn ta có : 
306
m 1
1
0
m m 2
1 2
0
m m m m n
1 2 n 1 n
0
cx
c
cx x
c
cx x x x
c−
⎧ ≈⎪⎪⎪ ≈⎪⎨⎪⋅ ⋅ ⋅⎪⎪ ⋅ ⋅ ⋅ ≈⎪⎩
              (6) 
Từ (6) ta có : 
m 1
1
0
m 2
2
1
m n
n
n 1
cx
c
cx
c
cx
c −
⎧ ≈⎪⎪⎪ ≈⎪⎨⎪⋅ ⋅ ⋅⎪⎪ ≈⎪⎩
                  (7) 
Ta suy ra trị tuyệt đối của các nghiệm là : 
1m1
0
2m2
1
nmn
n 1
cx
c
cx
c
cx
c −
⎧⎪ ≈⎪⎪⎪ ≈⎪⎨⎪⋅ ⋅ ⋅⎪⎪⎪ ≈⎪⎩
                (8) 
Sau đó ta thay ±|xi| vào (1) và tìm được nghiệm. Như vậy vấn đề là phải tìm 
được phương trình (3) dựa vào phương trình (1). Quá trình biến đổi (1) thành 
(3) được tiến hành như sau:  
Đầu tiên từ (1) ta xây dựng phương trình mới có nghiệm là  21x−  và quá 
trình này được gọi là bình phương nghiệm. Vì  P(x) có các nghiệm là xi nên nó 
có dạng: 
    P(x) = a0(x ‐ x1) (x ‐ x2)⋅⋅⋅( x ‐ xn) 
Do đó : 
    (‐1)nP(‐x) = a0(x + x1) (x + x2)⋅⋅⋅( x + xn) 
Ta suy ra : 
307
n 2 2 2 2 2 2 2
0 1 2 n( 1) P(x)P( x) a (x x )(x x ) (x x )− − = − − ⋅ ⋅ ⋅ −  
Thay x2 = ‐y ta có : 
2 2 2 2
1 0 1 2 nP (y) P(x)P( x) a (y x )(y x ) (y x )= − = + + ⋅ ⋅ ⋅ +  
Đa thức P1(y) có nghiệm là  2i 1y x= − . Đa thức này có dạng : 
(1) n (1) n 1 (1) n 2 (1)
1 0 1 2 nP (y) a y a y a y a
− −= + + + ⋅ ⋅ ⋅ +         (9) 
Do     P(x) = a0xn + a1xn‐1 + ⋅ ⋅ ⋅ + an = 0 
nên    P(‐x) = (‐1)n[a0xn ‐ a1xn‐1 + ⋅ ⋅ ⋅ + (‐1)nan] 
    (‐1)n P(‐x) = a0xn ‐ a1xn‐1 + ⋅ ⋅ ⋅ + (‐1)nan 
và            
n 2 2n 2 2n 2 2 2n 4
0 1 0 2 2 1 3 0 4
n 2
n
( 1) P(x)P( x) a x (a 2a a )x (a 2a a 2a a )x
( 1) a
− −− − = − − + − +
+ ⋅ ⋅ ⋅ + −  
Thay x2 = ‐y ta có : 
     2 n 2 n 1 2 n 2 n 21 0 1 0 2 2 1 3 0 4 nP (y) a y (a 2a a )y (a 2a a 2a a )y ( 1) a
− −= − − + − + + ⋅ ⋅ ⋅ + −       (10) 
So sánh (9) và (10) ta có công thức tính các hệ số của P1(y) từ các ak là : 
(1) 2
0 0
(1) 2
1 1 0 2
(1) 2
2 2 1 3 0 4
k
(1) 2 i
k k k i k i
i 1
(1) 2
n 1 n 1 n 2 n
(1) 2
n n
a a
a a 2a a
a a 2a a 2a a
a a 2 ( 1) a a
a a 2a a
a a
− +
=
− − −
⎧ =⎪ = −⎪⎪ = − +⎪⋅ ⋅ ⋅⎪⎪⎨ = + −⎪⎪⋅ ⋅ ⋅⎪⎪ = −⎪⎪ =⎩
∑             (11) 
Tiếp tục quá trình bình phương nghiệm đa thức P1(y) ta được P2(y) có 
nghiệm là  2 2k 1y ( x )= −  với các hệ số  (2)ka được tính theo  (1)ka  tương tự như (11) 
khi tính  (1)ka theo ak. Tiến hành bình phương nghiệm (s + 1)  lần ta có (3). Lúc 
đó các hệ số được xác định bằng: 
308
2(s 1) (s)
0 0
k2(s 1) (s) i (s) (s)
k k k i k i
i 1
2(s 1) (s) (s) (s)
n 1 n 1 n 2 n
2(s 1) (s)
n n
a a
a a 2 ( 1) a a
a a 2a a
a a
+
+
− +
=
+
− − −
+
⎧ ⎡ ⎤= ⎣ ⎦⎪⋅ ⋅ ⋅⎪⎪⎪ ⎡ ⎤= + − −⎣ ⎦⎪⎨⋅ ⋅ ⋅⎪⎪ ⎡ ⎤= −⎪ ⎣ ⎦⎪⎪ ⎡ ⎤= ⎣ ⎦⎩
∑
Vậy khi s đủ lớn ta có : 
s
(s)
2 m i
i i (s)
i 1
ax x
a −
− = − ≈  
Ta có m = 26 = 64. Dùng logarit ta tính ra các nghiệm : x1  = ‐4, x2 = 2, x3 = 1 
Ta xây dựng hàm  lobachevskygraeffe() để thực hiện thuật toán trên 
function y = lobachevskygraeffe(a, maxiter) 
% giai pt bang pp Lochevsky ‐ Graeffe 
% a la da thuc can tim nghiem 
% maxiter la so lan lap max 
c = a; 
n = length(a); 
m = 1; 
while m < maxiter 
    b = a; 
    for k = 2:n‐1; 
        s = 0; 
        i = 1; 
        while (((k ‐ i) >= 1)&((k + i) <= n))   
            s = s + (‐1)^i*b(k ‐ i)*b(k + i); 
            i = i + 1; 
        end 
        a(k) = b(k)^2 + 2*s; 
    end 
    a(1) = a(1)^2; 
    a(n) = a(n)^2; 
    j = 2^m; 
309
    for i = 1:n‐1 
        err = 1; 
        x(i) = a(i + 1)/a(i); 
        if x(i) == 1 
            x(i) = 1; 
        else 
            x(i) = exp((1/j)*log(x(i))); 
        end 
        err = abs(horner(c, x(i))); 
    end 
    if err < eps 
        break; 
    end 
    m = m + 1; 
end 
for i = 1:n‐1 
    if round(polyval(c, x(i))) ~= 0 
        x(i) = ‐x(i); 
    end 
end 
y = x; 
Để tìm nghiệm của đa thức P4(x) = x4 + 2x3 ‐ 25x2 ‐ 26x2 + 120 ta dùng chương 
trình ctlobachevskygraeffe.m: 
clc, clear all 
a = [1 2 ‐25 ‐26 120]; 
x = lobachevskygraeffe(a, 50) 
§20. PHƯƠNG PHÁP SCHRODER 
  Phương pháp lặp Schroder dùng để tìm nghiệm bội, có dạng tương tự 
như công thức lặp Newton: 
  + = − ′
k
k 1 k
k
mf(x )x x
f (x )
Trong đó m  là bội của nghiêm. Ban đầu  ta có  thể chưa biết m nên cần phải 
thử. Ta xây dựng hàm schroder() để thực hiện thuật toán trên: 
310
function [r, iter] = schroder(f1, df1, m, x0, tol) 
%Ham tim nghiem boi bang thuat toan Schroder 
iter = 0; 
d = feval(f1, x0)/feval(df1, x0); 
while abs(d) > tol 
    x1 = x0 ‐ m*d; 
    iter = iter + 1; 
    x0 = x1; 
    d = feval(f1, x0)/feval(df1, x0); 
end 
r = x0; 
Để giải phương trình  
  − − =x 2(e x) 0  
Ta dùng chương trình ctschroder.m với m = 2: 
clear all, clc 
[x, iter] = schroder(ʹf1ʹ, ʹdf1ʹ, 2, ‐2, 1e‐4) 
Trong đó: 
  function y = f1(x) 
y = (exp(‐x) ‐ x).^2; 
function y = df1(x) 
y = 2.0*(exp(‐x) ‐ x).*(‐exp(‐x) ‐ 1);