241
CHƯƠNG 5: CÁC PHƯƠNG TRÌNH PHI TUYẾN
§1.KHÁINIỆMCHUNG
Nếuphươngtrìnhđạisốhaysiêuviệtkháphứctạpthìítkhitìmđược
nghiệmđúng.Bởivậyviệctìmnghiệmgầnđúngvàướclượngsaisốlàrất
c
ầnthiết.
Taxétphươngtrình:
f(x)=0(1)
vớif(x)làhàmchotrướccủabiếnx.Chúngtacầntìmgiátrịgầnđúngcủa
nghiệmcủaphươngtrìnhnày.
Quátrìnhgi
ảithườngchialàmhaibước:bướcsơbộvàbướckiệntoàn
nghiệm.
Bước giải sơ bộ có 3 nhiệm vụ: vây nghiệm,tách nghiệmvàthuhẹp
khoảngchứanghiệm.
Vây
nghiệmlàtìmxemcácnghiệmcủaphươngtrìnhcóthểnằmtrên
nhữngđoạnnàocủatrụcx.Táchnghiệmlàtìmcáckhoảngchứanghiệmsao
cho trong mỗi khoảng chỉ cóđúng m
ột nghiệm. Thu hẹp khoảng chứa
nghiệmlà làmchokhoảngchứanghiệmcàngnhỏcàngtốt.Saubướcsơbộta
cókhoảngchứanghiệmđủnhỏ.Đểxácđịnhkhoảngchứanghiệmta
cóthể
dùngphươngphápđồthị.Ngoàiratacũngcóthểtìmnghiệmbằngphương
pháptìmtăngdần.Ýtưởngcủaphươngphápnàylànếyf
1(x).f2(x)<0thìcóít
nhấtmộtnghiệmcủaphươngtrìnhtrongđoạn[x
1,x2].Nếuđoạn[x1,x2]đủ
ệnthấykhoảng chứanghiệm,hàm trảvềgiátrịbiên củađoạn.
Nếukhôngcónghiệm,x1=x2=NaN.Tagọirootsearch()nhiềulầnđểphát
hiệnhếtcácđoạnchứanghiệm.
Vớivídụtìmkhoảngchứanghiệmcủahàm
f(x)=x
3
‐10x
2
+5tadùngchươngtrìnhctrootsearch.m
clearall,clc
f=inline(ʹx^3‐10*x^2+5ʹ);
[x1,x2]=rootsearch(f,2,10,.2)
Bướckiệntoànnghiệmtìmcácnghiệmgầnđúngtheoyêucầuđặtra.
Córấtnhiềuph
ươngphápxácđịnhnghiệmcủa(1).Sauđâychúngta
xéttừngphươngpháp.
§2.PHƯƠNGPHÁPLẶPĐƠN
Giảsửphươngtrình(1)đượcđưavềdạngtươngđương:
x=g(x)(2)
từgiátrịx
onàođógọilàgiátrịlặpđầutiêntalậpdãyxấpxỉbằngcôngthức:
x
n=g(xn‐1)(3)
vớin=1,2,
Hàmg(x)đượcgọilàhàmlặp.Nếudãyx
n→αkhin→∝thìtanóiphéplặp
(3)hộitụ.
%ra:x=nghiem
%err=
saiso|x(k)‐x(k‐1)|
%xx=cacgiatritrunggian
ifnargin<4
maxiter=100;
end
ifnargin<3
tolx=1e‐6;
end
xx(1)=x0;
fork=2:maxiter
xx(k)=feval(g,xx(k‐1));
err=abs(xx(k)‐xx(k‐1));
if
err<tolx
break;
end
x1
xo
x1
xo
244
end
x=xx(k);
ifk==maxiter
fprintf(ʹKhonghoitusau%dlanlap\nʹ,maxiter)
else
%a/b=biencuadoancantimnghiem
%tolx=saisomongmuon
%maxiterlanlapmax
%ra:x
=nghiem
%err=saiso
%xx=cacgiatritrunggian
y
x
a
b
ξ
b
1
245
tol=eps;
fa=feval(f,a);
fb=feval(f,b);
iffa*fb>0
error(ʹNghiemkhongotrongdoannayʹ);
end
fork=1:maxiter
xx(k)=(a+b)/2;
fx=feval(f,xx(k));
err=(b‐a)/2;
ifabs(fx)<tol|abs(err)<
tolx
1=a+h1
246
Trongđó
1
f(a)
(b a)
h
f(a) f(b)
=
−
−
−+
Tiếptheodùngcáchđóvớiđoạn[a,x
1]
hay[x
1,b]màhaiđầuhàmnhậngiátrịtrái
dấutađượcnghiệmgầnđúngx
2v.v.
Vềmặthìnhhọc,phươngphápnàycó
nghĩalàkẻdâycungcủađườngcongf(x)
quahaiđiểmA[a,f(a)]vàB[b,f(b)]haynóicáchkháclàtuyếntínhhoáhàm
f(x)trongđoạn[a,b].
Thật
vậyphươngtrìnhdâycungABcódạng:
f(a) f(b) af(b) bf(a)
yx
fa=feval(f,a);
fb=feval(f,b);
iffa*fb>0
error(ʹNghiemkhongotrongdoannay!ʹ);
end
fork=1:maxiter
xx(k)=(a*fb‐b*fa)/(fb‐fa);%pt.(1)
fx=feval(f,xx(k));
err=min(abs(xx(k)‐a),abs(b‐xx(k)));
a
b
x
y
x1
ξ
247
ifabs(fx)<tolfun|err<tolx
break;
elseiffx*fa>0
a=xx(k);
fa=fx;
else
b=xx(k);
fb=fx;
end
end
x=xx(k);
ifk==maxiter
0 f(x ) f (x )(x x ) O(x x )
++
′
=+ −+ −(2)
Giảsửrằngx
igầnvớixi+1,tacóthểbỏquasốhạngcuốitrong(2)vàcócông
thứcNewton‐Raphson:
i
i1 i
i
f(x )
xx
f(x)
+
=−
′
(3)
Nếux
i+1lànghiệmđúngcủaphươngtrìnhthìsaisốlàei=x‐xi.Khinghiệm
đượctínhtheo(3)thìsaisốlà:
248
2
i
i1 i
i
f(x)
ee
%ra:x=nghiem
%fx=f(x(last)),xx=cacgiatritrunggian
h=1e‐4;
h2=2*h;
tolf=eps;
ifnargin==4&isnumeric(df)
maxiter=tolx;
tolx=x0;
x0=df;
end
xx(1)=x0;
fx=feval(f,x0);
fork=1:maxiter
if~isnumeric(df)
dfdx=feval(df,xx(k));%daohamcuaham
else
dfdx=(feval(f,xx(k)+h)‐feval(f,xx(k)‐h))/h2;%daohamso
a
b
=xo
x1
y
x
249
end
dx=‐fx/dfdx;
xx(k+1)=xx(k)+dx;%pt.(3)
fx=feval(f,xx(k+1));
−
−
′
≈
−
(1)
vàtốnítthờigiantínhhơnkhidùngđạohàmgiảitíchhayđạohàmsố.
Bằngcáchxấpxỉ,biểuthức:
k
k1 k
k
f(x )
xx
f(x )
+
=−
′
trởthành:
−
+
−
⎡⎤
−
=− =−
⎢⎥
−
⎣⎦
Phéplặpcóthểkhônghộitụ(hìnhb).Tuynhiênkhihộitụ,nóhộirấtnhanh.
Taxâydựnghàm
secant()đểthựchiệnthuậttoántrên.
function[x,fx,xx]=secant(f,x0,x1,tolx,maxiter)
%giaiptf(x)=0bgphuongphapdaycung
%vao:f‐hamcantimnghiem
%x0,x1‐giatrikhoidongpheplap
%tolx‐saisomongmuon
%maxiter‐solanlapmax
%ra:x
=nghiem
%fx=f(x(last)),xx=cacgiatritrunggian
h=(x1‐x0)/100;
h2=2*h;
tolfun=eps;
xx(1)=x0;
fx=feval(f,x0);
fork=1:maxiter
ifk<=1
dfdx=(feval(f,xx(k)+h)‐feval(f,xx(k)‐h))/h2;
else
lap\nʹ,maxiter)
else
fprintf(ʹHoitusau%dlanlap\nʹ,k)
end
Đểtính nghiệm của hàm f(x) = x
3
‐10x
2
+ 5 ta dùng chương trình
ctsecant.m
clearall,clc
f=inline(ʹx.^3‐10*x.^2+5ʹ);
[x,ss,xx]=secant(f,0.5,1,1e‐4,50)
§7.PHƯƠNGPHÁPBRENT
Phương pháp Brent kêt hợp phương pháp chiađôi cung và phương
phápnộisuybậchaiđểtạo ra mộtphươngpháptìmnghiệm của phương
trìnhf(x) =0rấthiệuquả.Phương
phápnàydùngkhiđạohàmcủaf(x)khó
tínhhaykhôngthểtínhđược.Giảsửtacầntìmnghiêmtrong đoạn[x
1,x2].
Quátrìnhtìmnghiệmbắtđầubằngviệcchiađôiđoạn[x
1,x2]bằngđiểmx3.
−− −− −−
=++
−− −− −−
Choy=0tacó:
2312 3 1323 1 1231 2
122331
ffx(f f) ffx(f f) ffx(f f)
x
(f f )(f f )(f f )
−+ −+ −
=−
−−−
(1)
Độthayđổicủanghiệmlà:
31 2 2 3 1 212 3 123 1
33
122331
x (f f )( f f f ) f x (f f ) f x (f f )
xxx f
(f f )(f f )(f f )
−−++ −+ −
∆= − =
−−−
(2)
Taxâydựnghàm
brent()đểthựchiệnthuậttoántrên
fori=1:30
f3=feval(f,x3);
ifabs(f3)<tol
root=x3;
return
end
%xacdinhdoanchuanghiem.
if
f1*f3<0.0;
b=x3;
else
a=x3;
end
if(b‐a)<tol*max(abs(b),1.0)
root=0.5*(a+b);
return
end
%noisuybac2.
denom=(f2‐f1)*(f3‐f1)*(f2‐f3);
numer=x3*(f1‐f2)*(f2‐f3+f1)
+f2*x1*(f2‐f3)+
f1*x2*(f3‐f1);
ifdenom==0;
dx=b‐a;
else
dx=f3*numer/denom;
end
x=x3+dx;
%neulaprangoaidoan(a,b),dungcachchiadoicung
if(b‐x)*(x‐a)<0.0
Nhưvậ
y:
x
n+1=f(xn)(3)
x
n=f(xn‐1)(4)
Trừ(3)cho(4)vàápdụngđịnhlíLagrangechovếphảivớic∈[a,b]tacó:
x
n+1‐xn=f(xn)‐f(xn‐1)=(xn‐xn‐1)f’(c)(5)
Vìphéplặp(1)nên:
|x
n+1‐xn|≤q|xn‐xn‐1|(6)
Do(6)đúngvớimọinnênchon=1,2,3,...tacó:
|x
2‐x1|≤q|x1‐xo|
|x
3‐x2|≤q|x2‐x1|
255
...................
|x
n+1‐xn|≤q|xn‐xn‐1|
Điềunàycónghĩalàdãyx
i+1‐xi,mộtcáchgầnđúng,là mộtcấpsốnhân.Ta
cũngcoirằngdãyx
n‐yvớiylànghiệmđúngcủa(1),gầnđúngnhưmộtcấp
sốnhâncócôngsaiq.Nhưvậy:
n1
n
(10)
Thaygiátrịcủaqvừatínhở(10)vàobiểuthứccủaqởtrêntacó:
(
)
2
nn1
n
nn1n2
xx
yx
x2x x
+
++
−
=−
−+
(11)
Công thức (11)được gọi là công thức ngoại suy Adam. Nhưvậytheo(11)
trướchếttadùngphươngpháplặpđểtínhgiátrịgầnđúngx
n+2,xn+1,xncủa
nghiệmvàsauđótheo(11)tìmđượcnghiệmvớisaisốnhỏhơn.
Khi phương pháp lặpđược k ết hợp với phương pháp Aitken ta có
phươngphápSteffensen.Bắtđầulặptừx
0,haibướclặpđượcdùngđểtính:
x
1=f(x0) x2=f(x1)
vàsauđódùngthuậttoánAitkenđểtinhy
0theo(11).Đểtiếptụclặptacho
x
<maxiter)&(abs(f0‐x2)>.0000001))
count=count+1;
x1=feval(g,x0);
x2=feval(g,x1);
f0=x0‐(x1‐x0)*(x1‐x0)/(x2‐2.*x1+x0);
x0=f0;
end
y=f0;
Đểtìm nghiệmcủaphươngtrìnhx=(2‐e
x
+x
2
)/3=g(x)tadùngchươngtrình
ctaitstef.m
clearall,clc
f=inline(ʹ(2.‐exp(x)+x.^2)/3ʹ);
[x,y]=aitstef(f,0.,1e‐4,50)
§9.PHƯƠNGPHÁPMUELLER
Trongphươngphápdâycungkhitìmnghiệmtrongđoạn[a,b]taxấp
xỉ hàm bằng mộtđường thẳng. Tuy nhiênđểgiảm lượng tính toán vàđể
nghiệmhộitụnhanhhơntacóthểdùngphương
phápMuller.Nộidungcủa
phươngphápnàylàthayhàmtrongđoạn[a,b]bằngmộtđườngcongbậc2
màtahoàntoàncóthểtìmnghiệmchínhxáccủanó.
257
Gọicácđiểmđócóhoànhđộlầnlượtlàa=x2,b=x1vàtachọnthêm
2
22 22 2
v0(xx) a(0) b(0)cf
vh(xx) ah bh cf
vh(xx)ahbhcf
== ++=
== ++=
=− = − + =
Từđótacó:
0
1
2
101
2
1
201
fc
h
ahff
b
)1(h
f
)1(
f
f
a
=
−−
%giaiptf(x)=0
h1
h2
f(x)
x0,
f
0
x1,
f
1
x2,
f
2
258
%vao‐flahamcantimnghiem
%‐a,bladoanchuanghiem
%‐maxitersolanlapmax
%ra‐pxapxiMullercuaf
%‐ylagiatriy=f(p)
%‐errsaisothuccuanghiem.
%khoigana,b,x0vacacgiatrituongung
cuaham
x0=(a+b)/2.;
P=[x0ab];
Y=feval(f,P);
delta=1e‐4;
%tinhcachesocuaptbachai
fork=1:maxiter
h0=P(1)‐P(3);
P=R;
Y=feval(f,P);
end
%cactritinhlanmoi
P(3)=p;
Y(3)=feval(f,P(3));
y=Y(3);
%dieukien
dunglap
err=abs(z);
relerr=err/(abs(p)+delta);
if(err<delta)|(relerr<delta)|(abs(y)<eps)
break
end
end
Đểgiảiphươngtrìnhsin(x)‐0.5*x=0tadùngchươngtrình
ctmuller.m
clearall,clc
formatlong
f=inline(ʹsin(x)‐0.5*xʹ);
x=muller(f,1.8,2.2,50)
§10.PHƯƠNGPHÁPHALLEY
PhéplặpNewtontìmnghiệmcủahàmphươngtrìnhx=g(x)là:
′
f(x)
%vao:‐hamcantimnghiemf
%‐giatridaux0
%‐solanlapcucdaimaxiter
%ra‐nghiemx
%dungdaohamso
i=0;
h=1e‐4;
while(i
<maxiter)
f1=feval(f,x);
df1=(feval(f,x+h)‐feval(f,x‐h))/(2*h);
ddf1=(feval(f,x+h)‐2*feval(f,x)+feval(f,x‐h))/(h*h);
hx=x‐(f1/df1)*1./(1‐(f1*ddf1)/(2*(df1)^2))
x=hx;
i=i+1;
if(abs(f1)<eps)
break;
end
end
haydùnghàm
halley2()
functionx=halley2(f,x,maxiter)
%hamdungdetimnghiemcuaptf(x)=0
%vao:‐hamcantimnghiemf
%‐giatridaux
%‐solanlapcucdaimaxiter
%ra‐nghiemx
x=halley2(f,‐3,50)
§11.PHƯƠNGPHÁPCHEBYSHEV
Khitìmnghiệmcủaphươngtrìnhđạisốtuyếntínhhayphương trình
siêuviệtf(x)=0tacóthểdùngmộthàmcó4thôngsốđểxấpxỉhàmf(x)
y(x)=p
1+p2(x‐p3)
e
(1)
Cácthôngsốp
1 và p3tạosựchuyểndịchtheocáctrục;thôngsốpxácđịnh
biênđộvàecungcấpđộcongcủahàm.
Takhảosáthàmf(x)trênđoạn[a,b]trongđóf(a).f(b)<0,nghĩalàtrong
đoạn[a,b]
tồntạinghiệmcủaphươngtrìnhf(x) =0.Tacóthêmđiềukiện
fʹ(x).fʹʹ(x) ≠ 0 ∀x∈[a,b]. Gọi x
i∈[a,b]làlầnxấp xỉ thứ i của nghiệmthì
nghiệmlầnthứi+1theocôngthứcPopovskilà:
1
e
i1 i
2
fef.f
xx(e1)1 1
fe1f
+
⎡⎤
′′′
xx
f
+
−=−
′
vàđólàphéplặpNewton
Khie=0.5
2
2
iii
i1 i
33
ii
10.5f.f
f
f(x) f(x) f(x)
f
xx
ff(x)2f(x)
+
′′
+
⎛⎞
−
⎜⎟
′
′
×
df=
(feval(f,xx(k)+h)‐feval(f,xx(k)‐h))/h2;%daohamso
d2f=(feval(f,xx(k)+h)‐2*feval(f,xx(k))+feval(f,xx(k)‐h))/h^2;
dx=‐fx/df^3‐0.5*fx^2*d2f/df^3;
xx(k+1)=xx(k)+dx;
fx=feval(f,xx(k+1));
263
ifabs(fx)<tol|abs(dx)<tol
break;
end
end
x=xx(k+1);
Đểgiảiphươngtrìnhtadùngchươngtrình
ctchebyiter
clearall,clc
f=inline(ʹx.^3‐10*x.^2+5ʹ);
x=chebyiter(f,‐3,1e‐4)
§12.PHƯƠNGPHÁPNEWTONDÙNGCHOHỆPHITUYẾN
Phương pháp Newton có thể được tổng quát hoáđểgiải hệ phương
trìnhphituyếndạng:
1123 n
2123 n
n123 n
f (x ,x ,x , ,x ) 0
+
hay:
fʹ(x
i).∆x=‐f(xi)
với ∆x=x
i+1‐xi
Đốivớihệphươngtrình,côngthứclặplà:
J(X
i)∆X=‐F(Xi)
TrongđóJ(X
i)làtoántửJacobi.Nólàmộtmatrậnbậcn(n‐tươngứngvới
sốthànhphầntrongvectơX)códạng:
264
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⋅⋅⋅
∂
∂
∂
∂
∂
∂
=
n
n
3
n
2
n
1
n
n
2
3
2
2
2
1
2
n
1
3
1
2
1
PhươngphápNewtontuyếntínhhoáhệvànhưvậyvớimỗibướclặp
cầngiảimộthệphươngtrìnhtuyếntính(màbiếnlà∆X
i)xácđịnhbởicông
thứclặpchotớikhivectơX(x
1,x2,x3, ,xn)gầnvớinghiệm.
Taxâydựnghàm
new4sys()đểthựchiệnthuậttoánnày
function[P,iter,err]=new4sys(f,jf,P,max1)
%vao‐FlaheptluutrongM‐filef.m
%‐JFlamatranjacobiluutrongM‐filejf.m
%‐Pvectonghiembandau
%‐max1solanlapcucdai
%ra‐Plavetonghiem
%
‐itersolanlapthucte
%‐errsaiso
Y=f(P);
fork=1:max1
J=jf(P);
Q=P‐(J\Yʹ)ʹ;
Z=f(Q);
err=norm(Q‐P);
relerr=err/(norm(Q)+eps);
P=Q;
Y=Z;
iter=
k;
if(err<eps)|(relerr<eps)|(abs(Y)<eps)
break
functionf=f(p)
f=[(p(1)^2+p(1)*p(2)‐2),(2*p(1)*p(2)+p(2)^2‐3)];
Nộidungcủa
jf.m:
functionjf=jf(p)
jf=[(2*p(1)+p(2))p(1)
(2*p(1))(2*p(1)+2*p(2))];
Tacóthểdùnghàm
new4sys2()đểthựchiệnthuậttoán:
functionroot=new4sys2(func,x,maxiter)
%PhuongphapNewton‐Raphsondetimnghiem
%cuaheptfi(x1,x2, ,xn)=0,i=1,2, ,n.
%Cuphap:root=new4sys2(func,x,tol)
%vao:
%func=controham,trave[f1,f2, ,fn].
%x=vec
tobandau[x1,x2, ,xn].
%tol=saisomongmuon
%ra:
%root‐nghiemcuahe
ifsize(x,1)==1;
x=xʹ;
Copyright: Tài liệu đại học ©