Numerick\303\251 hled\303\241n\303\255 pevn\303\275ch bod\305\257 funkc\303\255Je t\305\231eba na\304\215\303\255st numerickou knihovnu, pop\305\231\303\255pad\304\233 otev\305\231\303\255t zkolabovanou sekci a na\304\215\303\255st pot\305\231ebn\303\251 p\305\231\303\255kazy odklepnut\303\255m.with(NumericalMethods);
# if not present run the next groupRoot:=proc(infun,inx,inmet,intol)
description "Solves f(x)=0 using classical numerical methods.",
"Root(f(x),xinit=,method=,tolerance=).",
"Optional options: stoppingcriterion=,output=,maxiterations=,maxnumber=,digits=",
"xinit=[<a>,<b>] for bisection,regfalsi,secant; otherwise xinit=<a>.",
"Option: method=bisection, =regfalsi, =secant, =newton, =illinois, =brent, =new, =fixedpoint.",
"Also method=[fixedpoint,<phi(x)>] for custom fixed point iterator.",
"Also method=[newton,<m>] for newton for multiplicity m: x_k-m*f(x_k)/f'(x_k).",
"Option: stoppingcriterion=absolute(default), =relative, =value, =bracket.",
"Option: output=information(default), =iterates, =value, =sequence.",
"Default maximal number of iterations is 50.",
"Default maximal allowed number in iterations is 1000000.",
"Default number of decimal digits shown is Digits.";
local lmet,leps,maxnum,maxiter,isseq,ldig,istwo,isbra,\134
# controls: method,tolerance,maxnumber,cutoff, output(sequence?iter?info?),digits,[a,b]?,bracket?
lfun,isreal,lnext,ltest,lprintiter,lprintiterC,lprintinf,lprintinfC,lprintbra,lprintintro,\134
# data: fixedpoint function,real or complex,iteration step procedure, stoptest procedure
# print procedures for iterates, information, bracketing, and introinfo (secant)
isfar,lseq,tmploc,kloc,iloc,\134
# loop control: loop end?,output,stopping test value,cycling vars
xloc,fxloc,yloc,fyloc,zloc,fzloc,xnew,fnew,tmplc,lspec,lspc,lflag,lmark;
# work registers (lspec,lspc,lflag,lmark for [newton,m],illinois,brent,new) (zloc fzloc for brent,new)
### check whether x assigned (messes up options)
if assigned(x) then
error("x has been assigned, can't use it as variable.") end if;
if nargs<4 then error("Some parameter(s) missing: f(x), xinit, method, tolerance.") end if;
### function(expression) input & check
if not(type(infun,algebraic)) or not(type(eval(infun,x=13.14),complex)) then
error("Expression with variable x expected as the first argument.") end if;
### method specification input & check & preparation
if not(type(inmet,equation)) or (convert(lhs(inmet),string)<>"method") then
error("The third argument must specify method=...") end if;
tmploc:=rhs(inmet);
if member(convert(tmploc,string),["bisection","regfalsi","secant","newton","fixedpoint","illinois","brent","new"],lmet) then
if lmet=1 then #bisection
lnext:=proc(xproc,yproc,fxp,fyp) return(evalf((xproc+yproc)/2)) end:
istwo:=true;isbra:=true
elif lmet=2 then #regfalsi
lnext:=proc(xproc,yproc,fxp,fyp)
if fyp=fxp then error("Division by zero, try different initial values.")
else return(evalf((xproc*fyp-yproc*fxp)/(fyp-fxp))) end if
end:
istwo:=true;isbra:=true
elif lmet=3 then #secant
lnext:=proc(xproc,yproc,fxp,fyp)
if fyp=fxp then error("Division by zero, try different initial values.")
else return(evalf((xproc*fyp-yproc*fxp)/(fyp-fxp))) end if
end:
istwo:=true;isbra:=false
elif lmet=4 then #newton
lfun:=infun;
lnext:=proc(fproc,xproc)
local dproc;
dproc:=eval(diff(fproc,x),x=xproc);
if dproc=0 then error("Division by zero, try different initial value.")
else return(evalf(xproc-eval(fproc,x=xproc)/dproc)) end if
end:
istwo:=false;isbra:=false
elif lmet=5 then #fixed point
lfun:='x'+infun;
lnext:=proc(fproc,xproc) return(evalf(eval(fproc,x=xproc))) end:
istwo:=false;isseq:=false
elif lmet=6 then #illinois
lnext:=proc(xproc,yproc,fxp,fyp)
if fyp=fxp then error("Division by zero, try different initial values.")
else return(evalf((xproc*fyp-yproc*fxp)/(fyp-fxp))) end if
end:
istwo:=true;isbra:=true
elif lmet=7 then #brent
lnext:=proc(xproc,yproc,fxp,fyp)
return(evalf((xproc*fyp-yproc*fxp)/(fyp-fxp)))
end:
istwo:=true;isbra:=true
elif lmet=8 then #new
lnext:=proc(xproc,yproc,tproc)
return(xproc+tproc*(yproc-xproc))
end:
istwo:=true;isbra:=true
end if
elif type(tmploc,list) and convert(tmploc[1],string)="fixedpoint" then #fixed point custom
lmet:=5;
lnext:=proc(fproc,xproc) return(evalf(eval(fproc,x=xproc))) end:
lfun:=tmploc[2];
if not(type(lfun,algebraic)) or not(type(eval(lfun,x=13.14),complex)) then
error("Custom fixed point expects an expression: method=[fixedpoint,<phi(x)>].")
end if;
if not(verify({fsolve(lfun=x,x)},{fsolve(infun=0,x)},('`subset`')(('float')(10^6)))) then
printf("The supplied iterative function seems to have a different fixed point than the root of the given function.\134n But you probably know what you are doing.\134n")
end if;
istwo:=false;isbra:=false
elif type(tmploc,list) and convert(tmploc[1],string)="newton" then #newton multiplicity m
lmet:=4;
lfun:=infun;
lspec:=tmploc[2];
if not(type(lspec,posint)) then
error("Modified Newton expects a positive integer: method=[newton,<m>].")
end if;
lnext:=proc(fproc,xproc)
local dproc;
dproc:=eval(diff(fproc,x),x=xproc);
if dproc=0 then error("Division by zero, try different initial value.")
else return(evalf(xproc-lspec*eval(fproc,x=xproc)/dproc)) end if
end:
istwo:=false;isbra:=false
else
error("Method not recognized, use method=bisection, =regfalsi, =secant, =newton, =illinois,\134
=brent, =new, =fixedpoint, =[fixedpoint,<phi(x)>], =[newton,<m>]")
end if;
### initial data input and check
isreal:=true; ### real numbers only, for printing format, necessary for some methods
if not(type(inx,equation)) or (convert(lhs(inx),string)<>"xinit")
then error("Second argument must specify initial value(s): xinit=<a> or xinit=[<a>,<b>].") end if;
tmploc:=rhs(inx);
if istwo then #two numbers expected
if type(tmploc,list(complex)) then
if isbra and not(type(tmploc,list(realcons))) then
error("For the chosen method, a list of two real initial values xinit=[<a>,<b>] is expected.")
end if;
xloc:=tmploc[1];
yloc:=tmploc[2];
fxloc:=evalf(eval(infun,x=xloc));
fyloc:=evalf(eval(infun,x=yloc));
if fxloc=0 then printf("Your initial guess is a root.\134n"); return(xloc) end if;
if fyloc=0 then printf("Your initial guess is a root.\134n"); return(yloc) end if;
if isbra then
if not(type(fxloc,realcons)) or not(type(fyloc,realcons))
or not(type(eval(infun,x=13.14),realcons)) then
error("Complex function values not allowed for the chosen method.")
end if
else
if not(type(xloc,realcons)) or not(type(fxloc,realcons))
or not(type(yloc,realcons)) or not(type(fyloc,realcons))
or not(type(eval(infun,x=13.14),realcons)) then isreal:=false
end if
end if;
if isbra and xloc>=yloc then
error("For the chosen method, xinit=[<a>,<b>] with a<b needed.")
end if;
if isbra and fxloc*fyloc>=0 then
error("For the chosen method, use xinit=[<a>,<b>] so that f(a), f(b) have opposite signs.")
end if
else
error("For the chosen method, a list of two initial values xinit=[<a>,<b>] is expected.")
end if
else
if type(tmploc,complex) then
xloc:=tmploc;
yloc:=0;
fxloc:=evalf(eval(infun,x=xloc));
fyloc:=0;
if fxloc=0 then printf("Your initial guess is a root.\134n"); return(xloc) end if;
if not(type(xloc,realcons)) or not(type(fxloc,realcons))
or not(type(eval(infun,x=13.14),realcons)) then isreal:=false
end if
else
error("For the chosen method, an initial number xinit=<a> is expected.")
end if
end if;
### tolerance input & check
if not(type(intol,equation)) or (convert(lhs(intol),string)<>"tolerance")
or not(type(evalf(rhs(intol)),positive)) then
error("Fourth argument must specify tolerance.") end if;
leps:=evalf(rhs(intol));
########################
### setting up defaults
isseq:=false;
maxiter:=50;
maxnum:=1000000;
ltest:=proc(fproc,xproc,lproc) ### default stopping condition absolute
return(abs(evalf(xproc-lproc)))
end;
ldig:=Digits;
############### printing routines
lprintbra:=proc(xproc,yproc)
printf(" a=% *.*f b=% *.*f ",ldig+4,ldig,xproc,ldig+4,ldig,yproc)
end:
lprintintro:=proc() printf("\134n") end:
lprintiter:=proc(xproc)
printf(" x=% *.*f",ldig+4,ldig,xproc)
end:
lprintiterC:=proc(xproc)
printf(" x=% *.*Zf",ldig+4,ldig,xproc)
end:
lprintinf:=proc(fproc,tproc)
printf(" f(x)=% *.*f test=% *.*f\134n",ldig+4,ldig,fproc,ldig+2,ldig,tproc)
end:
lprintinfC:=proc(fproc,tproc)
printf(" f(x)=% *.*Zf test=% *.*f\134n",ldig+4,ldig,fproc,ldig+2,ldig,tproc)
end:
if not(isreal) then
lprintiter:=lprintiterC;
lprintinf:=lprintinfC;
end if;
#######################################
### optional argument(s) processed here
if nargs>4 then
for kloc from 5 to nargs do
if not(type(args[kloc],equation)) then
error("The option %1 not recognized.
Optional options are stoppingcriterion=, output=, maxiterations=, maxnumber=, digits=.",args[kloc])
end if;
tmploc:=convert(lhs(args[kloc]),string);
if not(member(tmploc,["stoppingcriterion","output","maxiterations","maxnumber","digits"],'iloc')) then
error("The option %1 not recognized.
Optional options are stoppingcriterion=, output=, maxiterations=, maxnumber=, digits=.",args[kloc])
end if;
if iloc=5 then ######### digits
ldig:=rhs(args[kloc]);
if not(type(ldig,posint)) then
error("Number of decimal digits to show must be a positive integer.") end if
elif iloc=4 then ######### maxnumber
maxnum:=rhs(args[kloc]);
if not(type(maxnum,positive)) then
error("Maximal allowed number must be a positive number.") end if
elif iloc=3 then ######### maxiterations
maxiter:=rhs(args[kloc]);
if not(type(maxiter,posint)) then
error("Maximal number of iterations must be a positive integer.") end if
elif iloc=2 then ######### output
tmplc:=convert(rhs(args[kloc]),string);
if tmplc="sequence" then
lprintbra:=proc(xproc,yproc) end:
lprintintro:=proc() end:
lprintiter:=proc(xproc,yproc) end:
lprintiter:=proc(xproc) end:
lprintiterC:=proc(xproc) end:
lprintinf:=proc(fproc,tproc) end:
lprintinfC:=proc(fproc,tproc) end:
isseq:=true
elif tmplc="value" then
lprintbra:=proc(xproc,yproc) end:
lprintintro:=proc() end:
lprintiter:=proc(xproc,yproc) end:
lprintiter:=proc(xproc) end:
lprintiterC:=proc(xproc) end:
lprintinf:=proc(fproc,tproc) end:
lprintinfC:=proc(fproc,tproc) end:
isseq:=false
elif tmplc="information" then
isseq:=false
elif tmplc="iterates" then
lprintinf:=proc(fproc,tproc) printf("\134n") end:
lprintinfC:=proc(fproc,tproc) printf("\134n") end:
isseq:=false
else
error("Values for output are =information, =iterates, =value, =sequence.")
end if
else ######### stoppingcriterion
tmplc:=convert(rhs(args[kloc]),string);
if tmplc="absolute" then
ltest:=proc(fproc,xproc,lproc)
return(abs(evalf(xproc-lproc)))
end
elif tmplc="relative" then
ltest:=proc(fproc,xproc,lproc)
if xproc=0 then return(1) else return(abs(evalf((xproc-lproc)/xproc))) end if
end
elif tmplc="value" then
ltest:=proc(fproc,xproc,lproc)
return(abs(evalf(fproc)))
end
elif tmplc="bracket" then
if isbra then
ltest:=proc(fproc,xproc,lproc)
return(abs(evalf(yloc-xloc)))
end
else error("Stoppingcriterion=bracket only possible for bracketing methods.")
end if
else error("Values for stoppingcriterion are =absolute, =relative, =value, =bracket.")
end if
end if
end do
end if;
if ldig>Digits then Digits:=ldig end if;
##################################
isfar:=true;
xnew:=xloc;
if lmet<3 then ######### bisection or regfalsi
for kloc from 0 to maxiter by 1 while isfar do
printf("k=%02d ",kloc);
lprintbra(xloc,yloc);
tmplc:=xnew;
xnew:=lnext(xloc,yloc,fxloc,fyloc);
fnew:=evalf(eval(infun,x=xnew));
if not(type(fnew,realcons)) then
error("Complex function values appeared.") end if;
if isseq then lseq[kloc]:=xnew end if;
tmploc:=ltest(fnew,xnew,tmplc);
lprintiter(xnew);
lprintinf(fnew,tmploc);
if fnew=0 then isfar:=false;printf("Root (probably) found.") end if;
if tmploc<leps then isfar:=false end if;
if sign(fxloc)*sign(fnew)<0 then yloc:=xnew;fyloc:=fnew end if;
if sign(fnew)*sign(fyloc)<0 then xloc:=xnew;fxloc:=fnew end if;
if abs(fnew)>maxnum then maxnum:=0;isfar:=false end if
end do
elif lmet=3 then ######### secant
printf("k=%02d ",0);
lprintiter(xloc);lprintintro();
printf("k=%02d ",1);
lprintiter(yloc);lprintintro();
if isseq then lseq[0]:=xloc;lseq[1]:=yloc end if;
for kloc from 2 to maxiter by 1 while isfar do
xnew:=lnext(xloc,yloc,fxloc,fyloc);
xloc:=yloc;
yloc:=xnew;
fxloc:=fyloc;
fyloc:=evalf(eval(infun,x=yloc));
if not(type(yloc,realcons)) or not(type(fyloc,realcons)) then
isreal:=false;
lprintiter:=lprintiterC;
lprintinf:=lprintinfC;
end if;
if isseq then lseq[kloc]:=yloc end if;
tmploc:=ltest(fyloc,yloc,xloc);
printf("k=%02d ",kloc);
lprintiter(yloc);
lprintinf(fyloc,tmploc);
if fyloc=0 then isfar:=false;printf("Root (probably) found.") end if;
if tmploc<leps then isfar:=false end if;
if abs(yloc)>maxnum or abs(fyloc)>maxnum then maxnum:=0;isfar:=false end if
end do
elif lmet<6 then ######### newton, fixed point
printf("k=%02d ",0);
lprintiter(xloc);
lprintinf(fxloc,abs(fxloc));
for kloc from 1 to maxiter by 1 while isfar do
xnew:=lnext(lfun,xloc);
fnew:=evalf(eval(infun,x=xnew));
tmploc:=ltest(fnew,xnew,xloc);
xloc:=xnew;
if isseq then lseq[kloc]:=xnew end if;
if not(type(xnew,realcons)) or not(type(fnew,realcons)) then
isreal:=false;
lprintiter:=lprintiterC;
lprintinf:=lprintinfC;
end if;
printf("k=%02d ",kloc);
lprintiter(xloc);
lprintinf(fnew,tmploc);
if fnew=0 then
isfar:=false;
printf("Root (probably) found.")
end if;
if tmploc<leps then isfar:=false end if;
if abs(xnew)>maxnum or abs(fnew)>maxnum then maxnum:=0;isfar:=false end if
end do
elif lmet=6 then ######### illinois
lflag:=0;
lmark:=" ";
for kloc from 0 to maxiter by 1 while isfar do
printf("%c",lmark);printf("k=%02d ",kloc);
lprintbra(xloc,yloc);
tmplc:=xnew;
xnew:=lnext(xloc,yloc,fxloc,fyloc);
fnew:=evalf(eval(infun,x=xnew));
if not(type(fnew,realcons)) then
error("Complex function values appeared.") end if;
if isseq then lseq[kloc]:=xnew end if;
tmploc:=ltest(fnew,xnew,tmplc);
lprintiter(xnew);
lprintinf(fnew,tmploc);
if fnew=0 then isfar:=false;printf("Root (probably) found.") end if;
if tmploc<leps then isfar:=false end if;
lmark:=" ";
if sign(fxloc)*sign(fnew)<0 then # left end preserved
yloc:=xnew;fyloc:=fnew;
if lflag=-1 then fxloc:=fxloc/2;lmark:="M" end if;
lflag:=-1
end if;
if sign(fnew)*sign(fyloc)<0 then # right end preserved
xloc:=xnew;fxloc:=fnew;
if lflag=1 then fyloc:=fyloc/2;lmark:="M" end if;
lflag:=1
end if;
if abs(fnew)>maxnum then maxnum:=0;isfar:=false end if
end do
elif lmet=7 then ######### brent
if abs(fyloc)<abs(fxloc) then
zloc:=xloc;xloc:=yloc;yloc:=zloc;
fzloc:=fxloc;fxloc:=fyloc;fyloc:=fzloc;
end if;
zloc:=yloc;fzloc:=fyloc;
lflag:=true; # was bisection used?
for kloc from 0 to maxiter by 1 while isfar do
tmplc:=xnew;
if fxloc<>fzloc and fyloc<>fzloc then
xnew:=(xloc*fyloc*fzloc)/((fxloc-fyloc)*(fxloc-fzloc))
+(yloc*fxloc*fzloc)/((fyloc-fxloc)*(fyloc-fzloc))
+(zloc*fxloc*fyloc)/((fzloc-fxloc)*(fzloc-fyloc));
lmark:="Q"
else
xnew:=lnext(xloc,yloc,fxloc,fyloc);
lmark:="S"
end if;
lspc:=0;lspec:=(xloc+3*yloc)/4; # Choosing Bisection
if xloc<yloc and (xnew>lspec or xnew<=xloc) then lspc:=1 end if; # Is inside interval, not next to yloc?
if yloc<xloc and (xnew>=xloc or xnew<lspec) then lspc:=1 end if;
if lflag and abs(xnew-xloc)>=abs(xloc-zloc)/2 then lspc:=1 end if;
if not(lflag) and abs(xnew-xloc)>=abs(zloc-iloc)/2 then lspc:=1 end if;
if lflag and abs(xloc-zloc)<leps/2 then lspc:=1 end if;
if not(lflag) and abs(zloc-iloc)<leps/2 then lspc:=1 end if;
if lspc=1 then xnew:=(xloc+yloc)/2;lmark:="B";lflag:=true
else lflag:=false end if;
fnew:=evalf(eval(infun,x=xnew));
printf("%c",lmark);
printf("k=%02d ",kloc);
if xloc<yloc then lprintbra(xloc,yloc) else lprintbra(yloc,xloc) end if;
if not(type(fnew,realcons)) then
error("Complex function values appeared.") end if;
if isseq then lseq[kloc]:=xnew end if;
tmploc:=ltest(fnew,xnew,tmplc);
lprintiter(xnew);
lprintinf(fnew,tmploc);
if fnew=0 then isfar:=false;printf("Root (probably) found.") end if;
if tmploc<leps then isfar:=false end if;
iloc:=zloc;zloc:=xloc;fzloc:=fxloc;
if sign(fxloc)*sign(fnew)<0 then yloc:=xnew;fyloc:=fnew end if;
if sign(fyloc)*sign(fnew)<0 then xloc:=xnew;fxloc:=fnew end if;
if abs(fyloc)<abs(fxloc) then
lspc:=xloc;xloc:=yloc;yloc:=lspc;
lspc:=fxloc;fxloc:=fyloc;fyloc:=lspc;
end if;
if abs(fnew)>maxnum then maxnum:=0;isfar:=false end if;
end do
elif lmet=8 then ######### new (Chandrupatla)
zloc:=yloc;fzloc:=fyloc;
iloc:=0.5;isreal:=false;
lmark:="B";
for kloc from 0 to maxiter by 1 while isfar do
printf("%c",lmark);printf("k=%02d ",kloc);
if xloc<yloc then lprintbra(xloc,yloc) else lprintbra(yloc,xloc) end if;
tmplc:=xnew;
xnew:=lnext(xloc,yloc,iloc);
fnew:=evalf(eval(infun,x=xnew));
if not(type(fnew,realcons)) then
error("Complex function values appeared.") end if;
if isseq then lseq[kloc]:=xnew end if;
tmploc:=ltest(fnew,xnew,tmplc);
lprintiter(xnew);
lprintinf(fnew,tmploc);
if fnew=0 then isfar:=false;printf("Root (probably) found.") end if;
if tmploc<leps then isfar:=false end if;
if sign(fnew)*sign(fyloc)<0 then zloc:=xloc;fzloc:=fxloc end if;
if sign(fxloc)*sign(fnew)<0 then
zloc:=yloc;fzloc:=fyloc;yloc:=xloc;fyloc:=fxloc;
end if;
xloc:=xnew;fxloc:=fnew;
if abs(fnew)>maxnum then maxnum:=0;isfar:=false end if;
if abs(fxloc)<abs(fyloc) then lflag:=xloc else lflag:=yloc end if;
lflag:=(2*10^(1-Digits)+leps*abs(lflag))/abs(yloc-zloc);
lspc:=(xloc-yloc)/(zloc-yloc);
lspec:=(fxloc-fyloc)/(fzloc-fyloc);
if lspec^2<lspc and (1-lspec)^2<1-lspc then isreal:=true else isreal:=false end if;
if isreal then
iloc:=(fxloc*fzloc)/((fyloc-fxloc)*(fyloc-fzloc))\134
+((zloc-xloc)*fxloc*fyloc)/((yloc-xloc)*(fzloc-fxloc)*(fzloc-fyloc));
lmark:="Q"
else
iloc:=0.5;
lmark:="B"
end if;
iloc:=min(1-lflag,max(lflag,iloc))
end do
else
error("some problem with method name in execution block.")
end if;
if isfar then printf("Maximal number of iterations reached.") end if;
if maxnum=0 then printf("Numbers in iteration too large.") end if;
if isseq then
return(seq(lseq[iloc],iloc=1..kloc-1))
else
return(xnew)
end if
end:V dom\303\241c\303\255m \303\272kolu hled\303\241me ko\305\231en funkce x 2 - 3x + 1 p\305\231evodem na \303\272lohu pevn\303\251ho bodu.myf:=x^2-3*x+1;M\304\233li jste za \303\272kol naj\303\255t dva p\305\231evody a pro n\304\233 spo\304\215\303\255tat p\303\241r iterac\303\255. N\303\241sleduj\303\255c\303\255 k\303\263d v\303\241m umo\305\276n\303\255 zadat va\305\241i itera\304\215n\303\255 funkci jako myphi (te\304\217 tam je funkce ze standardn\303\255ho p\305\231evodu, kter\303\275 jste snad nepou\305\276ili). Vyzkou\305\241ejte, zda \304\215\303\255sla, kter\303\241 jste spo\304\215\303\255tali ru\304\215n\304\233, vych\303\241zej\303\255 i zde. M\304\233li jste tak\303\251 odhadnout, zda tato iterace bude \303\272sp\304\233\305\241n\303\241. Te\304\217 se uvid\303\255.myphi:=x^2-2*x+1:
Root(myf,xinit=1.,method=[fixedpoint,myphi],tolerance=0.0001);Ted si pohrajeme s relaxac\303\255.Uva\305\276ujeme \303\272lohu LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbW9HRiQ2LVEhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYjNiZGKy1GIzYmLUkmbWZyYWNHRiQ2KC1GIzYkLUkjbW5HRiQ2JFEiMkYnRi9GLy1GIzYmRistRiM2JC1JJW1zdXBHRiQ2JS1JI21pR0YkNiVRInhGJy8lJ2l0YWxpY0dGNy9GMFEnaXRhbGljRidGUS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGL0YrRi8vJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmZvLyUpYmV2ZWxsZWRHRjQtRiw2LVEiPUYnRi9GMi9GNkY0RjhGOkY8Rj5GQC9GQ1EsMC4yNzc3Nzc4ZW1GJy9GRkZgcEZmbkYvRitGL0YrRi8=. P\305\231\303\255slu\305\241nou itera\304\215n\303\255 funkci myphi jsme zadali n\303\255\305\276e a p\305\231idali upravenou funkci odpov\303\255daj\303\255c\303\255 ko\305\231enov\303\251 formulaci t\303\251to \303\272lohy.myphi:=2/x^2;
myf:=2/x^2-x;Jako obvykle hodl\303\241me za\304\215\303\255t v x 0 = 1. V k\303\263du n\303\255\305\276e je p\305\231ipravena relaxovan\303\241 iterace, nejprve s hodnotou LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbW9HRiQ2LVEhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYjNiZGKy1GIzYmLUkjbWlHRiQ2JVEnJiM5NTU7RicvJSdpdGFsaWNHRjRGLy1GLDYtUSI9RidGL0YyL0Y2RjRGOEY6RjxGPkZAL0ZDUSwwLjI3Nzc3NzhlbUYnL0ZGRlctSSNtbkdGJDYkUSIxRidGL0YvRitGL0YrRi8=, kdy vlastn\304\233 o relaxaci ani nejde.Spou\305\241t\304\233jte relaxovanou iteraci s parametrem LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbW9HRiQ2LVEhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYjNiQtSSNtaUdGJDYlUScmIzk1NTtGJy8lJ2l0YWxpY0dGNEYvRi9GK0Yv po\304\215\303\255naje hodnotou 1 a po 0.1 dosk\303\241kejte a\305\276 po 0.1, p\305\231itom sledujte, co to d\304\233l\303\241 s konvergenc\303\255 a jak rychle (za kolik krok\305\257) algoritmus dojede k c\303\255li. Je dobr\303\251 si p\305\231ipravit tabulku, kam budete pro testovan\303\251 lambdy ud\303\241vat po\304\215et krok\305\257 iterace, pop\305\231\303\255pad\304\233 N pro divergenci.lamb:=0.3:
Root(myf,xinit=1.,method=[fixedpoint,lamb*myphi+(1-lamb)*x],tolerance=0.0001);Nyn\303\255 sestavte vzorec pro relaxovanou itera\304\215n\303\255 funkci fi lambda a pomoc\303\255 jeho derivace najd\304\233te optim\303\241ln\303\255 hodnotu parametru lambda pro okol\303\255 bodu x 0 = 1. Souhlas\303\255 to s v\303\275sledky va\305\241ich experiment\305\257?Prost\303\241 iterace diverguje, co\305\276 plat\303\255 i pro hodnoty lambdy bl\303\255zk\303\251 jedni\304\215ce. Pak ale nastoup\303\255 konvergence a nejl\303\251pe to vypad\303\241, kdy\305\276 je lambda okolo 0.3.Optim\303\241ln\303\255 lambda: M\304\233la by vyj\303\255t itera\304\215n\303\255 funkce LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbW9HRiQ2LVEhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYjNiZGKy1GIzYoRistRiM2Ji1JI21pR0YkNiVRJyYjOTU1O0YnLyUnaXRhbGljR0Y0Ri8tRiw2LVExJkludmlzaWJsZVRpbWVzO0YnRi9GMi9GNkY0RjhGOkY8Rj5GQEZCL0ZGRkQtSSZtZnJhY0dGJDYoLUYjNiQtSSNtbkdGJDYkUSIyRidGL0YvLUYjNiZGKy1GIzYkLUklbXN1cEdGJDYlLUZPNiVRInhGJy9GU0Y3L0YwUSdpdGFsaWNGJ0Zobi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGL0YrRi8vJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmFwLyUpYmV2ZWxsZWRHRjRGLy1GLDYtUSIrRidGL0YyRldGOEY6RjxGPkZAL0ZDUSwwLjIyMjIyMjJlbUYnL0ZGRmpwLUYjNiYtSShtZmVuY2VkR0YkNiQtRiM2Ji1GaW42JEZecEYvLUYsNi1RKCZtaW51cztGJ0YvRjJGV0Y4RjpGPEY+RkBGaXBGW3FGTkYvRi9GVEZjb0YvRitGL0YrRi9GK0Yv.Polo\305\276\303\255me jej\303\255 derivaci LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbW9HRiQ2LVEhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYjNiZGKy1GIzYoLUYsNi1RKiZ1bWludXMwO0YnRi9GMi9GNkY0RjhGOkY8Rj5GQC9GQ1EsMC4yMjIyMjIyZW1GJy9GRkZRLUYjNiYtSSNtaUdGJDYlUScmIzk1NTtGJy8lJ2l0YWxpY0dGNEYvLUYsNi1RMSZJbnZpc2libGVUaW1lcztGJ0YvRjJGT0Y4RjpGPEY+RkBGQi9GRkZELUkmbWZyYWNHRiQ2KC1GIzYkLUkjbW5HRiQ2JFEiNEYnRi9GLy1GIzYmRistRiM2JC1JJW1zdXBHRiQ2JS1GVjYlUSJ4RicvRlpGNy9GMFEnaXRhbGljRictRl9vNiRRIjNGJ0YvLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0YvRitGLy8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGanAvJSliZXZlbGxlZEdGNEYvLUYsNi1RIitGJ0YvRjJGT0Y4RjpGPEY+RkBGUEZSLUYjNiYtRl9vNiRGZ3BGLy1GLDYtUSgmbWludXM7RidGL0YyRk9GOEY6RjxGPkZARlBGUkZVRi9GK0YvRitGL0YrRi8= rovnou nule, polo\305\276\303\255me x rovno jedn\303\251 a najdeme LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbW9HRiQ2LVEhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYjNiZGKy1GIzYmLUkjbWlHRiQ2JVEnJiM5NTU7RicvJSdpdGFsaWNHRjRGLy1GLDYtUSI9RidGL0YyL0Y2RjRGOEY6RjxGPkZAL0ZDUSwwLjI3Nzc3NzhlbUYnL0ZGRlctSSZtZnJhY0dGJDYoLUYjNiQtSSNtbkdGJDYkUSIxRidGL0YvLUYjNiQtRmluNiRRIjVGJ0YvRi8vJS5saW5ldGhpY2tuZXNzR0Zbby8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zlby8lKWJldmVsbGVkR0Y0Ri9GK0YvRitGLw== = 0.2, co\305\276 zhruba odpov\303\255d\303\241 numerick\303\275m experiment\305\257m.Pokud bychom zkou\305\241eli d\304\233lat derivaci nulovou v bod\304\233 1.26, tedy u \305\231e\305\241en\303\255, dostali bychom optim\303\241ln\303\255 lambdu, kter\303\241 se velmi dob\305\231e shoduje s experimentem.