#                                                   c2
# Let us use the following precision of floating point numbers.
> Digits:=10:
# 


# Define the number of iterations for bisection and Newton's Newton's modified 
# method. The results on each iterations are kept in the vectors x and y.
> n:=10:
> x:=array(0..n): y:=array(0..n):z:=array(0..n):
# 


# This is the implementation of the Newton's method. The parameter "a" is the initial 
# approximation, "n" is the number of iterations, "sol" is the solution obtained after n 
# iterations, "f" is the function whose roots we are looking for and  "trace" is the 
# vector of the results at every iteration.
> newton:=proc(a,n,sol,f,trace)
>   local i,t:
>   der:=D(f):
>   trace[0]:=evalf(a):
>   t:=evalf(a):
>   for i from 1 to n do
>       t:=t-evalf(f(t))/evalf(der(t)):
>       trace[i]:=t:
>  od:
> sol:=t:
> end:
> 
# 
Warning, `der` is implicitly declared local



# This is the implementation of the bisection method. The parameters "a" and "b" are 
# the boundaries of the interval in which we seek the solution, "n" is the number of 
# iterations, "sol" is the solution obtained after n iterations, "f" is the function whose 
# roots we are looking for and  "trace" is the vector of the results at every iteration.
> bisection:=proc(a,b,n,solu,f, trace)
>    local c,i,u,v,w,a1,b1:
>    a1:=evalf(a):  b1:=evalf(b):
>    u:=f(a1): v:=f(b1):
>    if sign(u)=sign(v) then
>    lprint(`the equation may not have a solution in the interval [`,a,b,`]`):
>    else
>      for i from 1 to n do
>       c:=(a1+b1)/2: w:=f(c): trace[i]:=c:
>       if sign(u)<>sign(w) then b1:=evalf(c):  v:=w: else a1:=evalf(c):  u:=w: fi:
>     od:
>     solu:=(a1+b1)/2:
>    fi:
> end:
> 


# Telo nasledujucej procedury newtonmod obsahuje implementaciu Newtonovej 
# metody.
# Upravte proceduru newtonmod tak aby sa z nej stala implementacia 
# modifikovanej Newtonovej metody tj. 
#                                     x[n+1]:=x[n]-p(x[n])/p'(x[0])
# 
> newtonmod:=proc(a,n,solm,f,trace)
>   local i,t:
>   der:=D(f):
>   trace[0]:=evalf(a):
>   t:=evalf(a):
>   for i from 1 to n do
>       t:=t-evalf(f(t))/evalf(der(t)):
>       trace[i]:=t:
>  od:
> solm:=t:
> end:
Warning, `der` is implicitly declared local



# Najdite korene nasledujucich rovnic s presnostou 0.00001 pomocou 
#                                                      Newtonovej metoy,
#                                                      modifikovanej  Newtonovej metdoy,
#                                                      a metody bisekcie.
# 
# Vysledky skontrolujte pomocou MAPLE-prikazu fsolve.
# 
# Pouzite nasledujuci postup:
> 
# Digits:=?:
# n:=?:
# x:=array(0..n): y:=array(0..n):z:=array(0..n):
# p:=t->????????????:
# plot(p(t),t=?..?);
# bisection(?,?,n,'solb','p','y'):
# newton(?,n,'soln','p','x'):
# newtonmod(?,n,'solm','p','z'):
# lprint(`Newton=`,soln,`Newtonmod=`,solm,` Biseccion=`,solb):
# lprint(` Iteration     Newton           Newtonmod          Bisection`):
# for i from 1 to n do
# printf(`%2i      %15.13f     %15.13f    %15.13f`,i,x[i],z[i],y[i]):lprint():
# od;
# fsolve(p(t),t);
> 
# Korene zapiste do tabulky:
# Priklad                                         Riesenie
# 1
# 2
# 3
# 4
# 5
# 6
# 7
# 8
> 


# 1.
# 2.
# 3.
# 4.
# 5.
# 6.
# 7.
# 8.
# 9.
# 10.
# 11.
# Najdite horeuvedenym postupom korene Wilkisonovho polynomu
> w20:=t->(t-1.)*(t-2.)*(t-3.)*(t-4.)*(t-5.)*(t-6.)*(t-7.)*(t-8.)*(t-9.)*(t-10.)*(t-11.)*(t-12.)*(t-13.)
> *(t-14.)*(t-15.)*(t-16.)*(t-17.)*(t-18.)*(t-19.)*(t-20.);
# 12. A teraz korene perturbovaneho Wilkinsonv polynomu
> w20p:=t->w20(t)-2^(-20)*t^19;
# 13. Najdite VSETKY korene polynomov z Prikladu 11 a 12 pomocou 
# MAPLE-prikazu fsolve
> fsolve(w20(t),t);
> 
> vys1:=fsolve(w20(t),t):
> vys2:=fsolve(w20p(t),t):
> vys3:=fsolve(w20p(t),t,complex):
> 'vys1'=vys1;
> 'vys2'=vys2;
> 'vys3'=vys3;
# 14.                           Nasobny koren
# Hladajte newtonovou metodou korene rovnice 
# Aka je rychlost konvergencie k nasobnemu korenu t=1. ?
# Modifikujte newtonovu metodu tak aby ste dostali kvadraticku konvergenciu.
> newton2:=proc(a,n,sol,f,trace)
>   local i,t:
>   der:=D(f):
>   trace[0]:=evalf(a):
>   t:=evalf(a):
>   for i from 1 to n do
# TU TREBA NIECO ZMENIT 
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
>   t:=t-evalf(f(t))/evalf(der(t)):
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>       trace[i]:=t:
>  od:
> sol:=t:
> end:
> 
# Pouzite nasledujuci postup
> Digits:=?
> n:=?:
> x:=array(0..?): y:=array(0..?):z:=array(0..?):w:=array(0..?):
> p:=t->??????????
> newton(?,n,'soln','p','x'):
> newtonmod(?,n,'solm','p','z'):
> newton2(?,n,'soln','p','w'):
> lprint(`Newton=`,soln,`Newtonmod=`,solm,` Newton2=`,solb):
> lprint(` Iteration     Newton           Newtonmod          Newton2`):
> for i from 1 to n do
> printf(`%2i      %15.13f     %15.13f    %15.13f`,i,x[i],z[i],y[i]):lprint():
> od;
> fsolve(p(t),t);


# 14. Este zopar 'divokych' funkcii (pouzite Newtona a bisekciu)
#   na intervale <-1.15>
# na intervale <1,1.57>
# na intervale <-2,2>
> t^3-2*t+2+0.1*sin(1000*t);
# 
# na intervale <-2,3>
> t^3-2*t+0.5*sin(1/t);
> 15. Skuste najst priklad kde je bisekcia rychlejsia ako Newton (skuste nasobne 
> korene).