In [1]:
# Numerické riešenie nelineárnych rovníc
# Metóda polovičného delenia intervalu (bisekcie)
# Newton Raphsonova metóda (regula falsi, metóda dotyčníc)
#
# Programové kódy: SageMath v.9.x(Python 3.0)
# Prerekvizity: Diferenciálny počet, polynómy, Python, LaTex.
#
# Použitá literatúra:
# [1] R. E. Mezei, An introduction to SAGE programming, John Wiley,and Sons,2016
# ISBN: 978-1-119-12280-7.
#
# [2] Gregory V.Bard. Sage for Undergraduates:Second Edition,Compatible with Python3.
# American Mathematical Society,2022.ISBN978-1-4704-6155-3.
# URL https://bookstore.ams.org/mbk-143/ - Sage pre pokročilých
#
# [3] Paul, Zimmermann and Alexandre, Casamayou and Nathann, Cohen and Guillaume Con-
# nan and Thierry Dumont and Laurent Fousse and François Maltey and Matthias Meulien
# and Marc Mezzarobba and Clément Pernet and Nicolas M. Thiéry and Erik Bray and
# John Cremona and Marcelo Forets and Alexandru Ghitza and Hugh Thomas:
# Computational Mathematics with SageMath, 2018, ISBN978-1-611975-45-1.
# https://my.siam.org/Store/Product/viewproduct/?ProductId=30174951,
# https://www.sagemath.org/sagebook/english.html
#
# E-kurz: P. Szabó: Numerická matematika, 2022
# Spúšťanie programov:
# 1) Otvoriť lokalitu https://sagecell.sagemath.org/
# 2) Kopírovať obsah príkladu In[n] do SageMath boxu
# 3) Kliknutie na "Evaluate"
# Príloha publikácie:
# [4] Peter Szabó a Eva Baranová: Integrálny počet a numerická matematika, 2022, TUKE
# ISBN:978-80-553-4185-9
# Modifikované: 10.12.2022
In [2]:
# Bisekcia - jednoduchá aplikácia algoritmu
# Je nutné zistiť interval na ktorom sú predpoklady metódy splnené
@interact
def BisectionInteract(
f = input_box(default = x^2-7, label = "f(x)" ),
a = input_box(default = 0 ),
b = input_box(default = 7 ),
tolerance = input_box(default= 0.00001) ):
x0 = (a+b)/2
while( f(x0)!=0 and (b-a) >= tolerance):
if(f(a)*f(x0)<0):
b = x0
else: #(f(b)*f(x0)<0)
a = x0
x0 = (a+b)/2
print ("approximate solution = ", x0.n())
In [3]:
# Bisekcia I - porovnanie výsledkov
show(sqrt(7))
sqrt(7).n()
# je potrebné zadať väčšiu presnosť, toleranciu
Out[3]:
In [4]:
# vypočítať koreň funkcie f(x)=cos(x) na intervale [0,pi]
# vypočítať hodnotu pi/2
In [5]:
# Bisekcia - test podmienky (f(b)*f(x0)<0)
# zadajte a=3
@interact
def BisectionInteract(
f = input_box(default = x^2-7, label = "f(x)" ),
a = input_box(default = 0 ),
b = input_box(default = 7 ),
tolerance = input_box(default= 0.00001),
auto_update=False ):
f(x)=f
#validation
if f(a)*f(b)>=0:
print ("Error: f(a)*f(b) should be negative!")
else:
x0 = (a+b)/2
while( f(x0)!=0 and (b-a)>= tolerance):
if(f(a)*f(x0)<0):
b = x0
else: #(f(b)*f(x0)<0)
a = x0
x0 = (a+b)/2
print ("approximate solution = ", x0.n())
In [6]:
# Bisekcia - krokovanie algoritmu, medzivýsledky, intervaly [ai,bi]
@interact
def BisectionInteract(
f = input_box(default = x^2-7, label = "f(x)" ),
a = input_box(default = 0 ),
b = input_box(default = 7 ),
tolerance = input_box(default= 0.00001),
auto_update=False ):
f(x)=f
#validation
if f(a)*f(b)>=0:
print ("Error: f(a)*f(b) should be negative!")
else:
x0 = (a+b)/2
while( f(x0)!=0 and (b-a)>= tolerance):
print ("on [", round(a, 4) , "," ,round(b, 4) ,"]")
print ("--- x0 = ",round( x0, 4), ", f(x0)=", round( f(x0), 7))
if(f(a)*f(x0)<0):
b = x0
else: #(f(b)*f(x0)<0)
a = x0
x0 = (a+b)/2
print ("approximate solution = ", x0.n())
In [7]:
# Bisekcia - krokovanie algoritmu, medzivýsledky, grafické výstupy
@interact
def BisectionInteract(
f = input_box(default = x^2-7, label = "f(x)" ),
a = input_box(default = 0 ),
b = input_box(default = 7 ),
tolerance = input_box(default= 0.00001),
auto_update=False ):
#save these values for plotting
A = a
B = b
#suppresses a warning message
f(x) = f
#validation
if f(a)*f(b)>=0:
print ("Error: f(a)*f(b) should be negative!")
else:
x0 = (a+b)/2
#this list will contain successive approximations
myList = [(x0,0)]
while( f(x0)!=0 and (b-a)>= tolerance):
if(f(a)*f(x0)<0):
b = x0
else: #(f(b)*f(x0)<0)
a = x0
x0 = (a+b)/2
myList = myList + [(x0,0)]
#plot the list of points
p = list_plot(myList, color = "red", size = 50)
#plot the function
q = plot(f(x), x, A, B)
(p+q).show()
In [8]:
# Bisekcia - krokovanie algoritmu, medzivýsledky, grafické výstupy
# Vstupy:
f(x) = x^2-7
a=0
b=7
tolerance= 0.00001
#save these values for plotting
A = a
B = b
#validation
if f(a)*f(b)>=0:
print ("Error: f(a)*f(b) should be negative!")
else:
x0 = (a+b)/2
#this list will contain successive approximations
myList = [(x0,0)]
while( f(x0)!=0 and (b-a)>= tolerance):
if(f(a)*f(x0)<0):
b = x0
else: #(f(b)*f(x0)<0)
a = x0
x0 = (a+b)/2
myList = myList + [(x0,0)]
#plot the list of points
p = list_plot(myList, color = "red", size = 50)
#plot the function
q = plot(f(x), x, A, B)
(p+q).show()
In [9]:
# Newtonova Raphsonova metóda na riešenie rovníc - algoritmus
@interact
def NewtonRaphsonInteract(
f = input_box(default = x^2-7, label = "f(x)" ),
x1 = input_box(default = 3 ),
max_steps = input_box(default = 100 ),
tolerance = input_box(default= 0.00001) ):
#suppresses a warning message
f(x) = f
#we start with the initial value x1
xn = x1
#we need the derivative of f with respect to x
fprime(x) = derivative(f(x),x)
steps = 1
#keep computing successive approximations as long as
# we did not reach max_steps
# we do not have a reasonable solution and
# the derivative is not too close to zero
while( steps<=max_steps and \
abs(f(xn))>tolerance and \
abs(fprime(xn))>tolerance ):
xn -= f(xn)/fprime(xn)
steps += 1
if(abs(fprime(xn))<=tolerance ):
print ("error!!! derivative value is too close to zero", xn.n())
elif (steps>=max_steps):
print ("I gave up ... too many steps!")
else:
print ("approximate solution = ", xn.n())
In [10]:
# presný výpočet
sqrt(7).n(100)
Out[10]:
In [11]:
# x1=0
# metóda nie je konvergentná
In [12]:
# Newtonova Raphsonova metóda na riešenie rovníc - algoritmus a graf riešenia
@interact
def NewtonRaphsonInteract2(
f = input_box(default = x^2-7, label = "f(x)" ),
x1 = input_box(default = 3 ),
max_steps = input_box(default = 100 ),
tolerance = input_box(default= 0.00001),
auto_update=False ):
#suppresses a warning message
f(x) = f
#we start with x1
xn = x1
#this will collect the list of all successive approximations
myList = [(xn,0)]
#we need the derivative of f with respect to x
fprime(x) = derivative(f(x),x)
steps = 1
#keep computing successive approximations as long as
# we did not reach max_steps
# we do not have a reasonable solution and
# the derivative is not too close to zero
while( steps<=max_steps and \
abs(f(xn))>tolerance and \
abs(fprime(xn))>tolerance ):
xn -= f(xn)/fprime(xn)
myList += [(xn,0)]
steps += 1
#output the results
if(abs(fprime(xn))<= tolerance ):
print ("error, derivative value close to zero", xn.n())
elif (steps>=max_steps):
print ("I gave up ... too many steps!")
else:
print ("approximate solution = ", xn.n())
#make the plot either way - so it can be analyzed
# - plot the list of points
p = list_plot(myList, color = "red", size = 50)
# - plot the function
a = min([myList[i][0] for i in [0,1,..,len(myList)-1]])
b = max([myList[i][0] for i in [0,1,..,len(myList)-1]])
q = plot(f(x), x, a-5, b+5)
(p+q).show()
In [13]:
# Newtonova Raphsonova metóda na riešenie rovníc - iteračné kroky algoritmu
@interact
def NewtonRaphsonInteract3(
f = input_box(default = x^2-7, label = "f(x)" ),
x1 = input_box(default = 5 ),
iterations = slider(vmin=1, vmax=10,
step_size=1,
default=1),
tolerance = input_box(default= 0.00001) ):
#suppresses a warning message
f(x) = f
#we start with x1
xn = x1
#this will collect the list of all successive approximations
myList = [(xn,0)]
#we need the derivative of f with respect to x
fprime(x) = derivative(f(x),x)
steps = 1
#keep computing successive approximations as long as
# we did not reach max_steps
# we do not have a reasonable solution and
# the derivative is not too close to zero
while( steps<=iterations and abs(f(xn))>tolerance and abs(fprime(xn))>tolerance ):
xn -= f(xn)/fprime(xn)
myList += [(xn,0)]
steps += 1
#output the results
if(abs(fprime(xn))<= tolerance ):
print ("error, derivative value close to zero", xn.n())
else:
print ("approximate solution = ", xn.n())
#make the plot either way - so it can be analyzed
# - plot the list of points
p = list_plot(myList, color = "red", size = 50)
# - plot the function
a = min([myList[i][0] for i in [0,1,..,len(myList)-1]])
b = max([myList[i][0] for i in [0,1,..,len(myList)-1]])
q = plot(f(x), x, a-5, b+5)
(p+q).show()
In [14]:
# Newtonova Raphsonova metóda na riešenie rovníc - iteračné kroky algoritmu II
@interact
def NewtonRaphsonInteract4(
f = input_box(default = x^2-7, label = "f(x)" ),
x1 = input_box(default = 30 ),
iterations = slider(vmin=1, vmax=10,
step_size=1,
default=1),
tolerance = input_box(default= 0.00001),
auto_update=False ):
#suppresses a warning message
f(x) = f
#we start with x1
xn = x1
#this will collect the list of all successive approximations
myList = [(xn,0)]
#we need the derivative of f with respect to x
fprime(x) = derivative(f(x),x)
steps = 1
#keep computing successive approximations as long as
# we did not reach max_steps
# we do not have a reasonable solution and
# the derivative is not too close to zero
while( steps<=iterations and \
abs(f(xn))>tolerance and \
abs(fprime(xn))>tolerance ):
xn -= f(xn)/fprime(xn)
myList += [(xn,0)]
steps += 1
#output the results
if(abs(fprime(xn))<= tolerance ):
print ("error, derivative value close to zero", xn.n())
else:
print ("approximate solution = ", xn.n())
#make the plot either way - so it can be analyzed
# - plot the list of points
p = list_plot(myList, color = "red", size = 50)
# - to put the text, we need to find how big is the graph
y1 = min([f(myList[i][0]) for i in [0,1,..,len(myList)-1]])
y2 = max([f(myList[i][0]) for i in [0,1,..,len(myList)-1]])
posT = y1-0.10*(y2-y1)
t = text(1, (myList[0][0], posT), fontsize=20, color = "black")
for i in range(1, len(myList)):
t = t+ text(i+1, (myList[i][0],posT), fontsize=20, color = "black")
# - plot the function
a = min([myList[i][0] for i in [0,1,..,len(myList)-1]])
b = max([myList[i][0] for i in [0,1,..,len(myList)-1]])
q = plot(f(x), x, a-5, b+5)
(p+q+t).show()
In [15]:
# zadajte funkciu f(x)=cos(x),a x1=25
In [16]:
@interact
def NewtonRaphsonInteract5(
f = input_box(default = x^2-7, label = "f(x)" ),
x1 = input_box(default = 30 ),
iterations = slider(vmin=1, vmax=10,
step_size=1,
default=1),
tolerance = input_box(default= 0.00001),
auto_update=False ):
#suppresses a warning message
f(x) = f
#we start with x1
xn = x1
#this will collect the list of all successive approximations
myList = [(xn,0)]
#we need the derivative of f with respect to x
fprime(x) = derivative(f(x),x)
steps = 1
#keep computing successive approximations as long as
# we did not reach max_steps
# we do not have a reasonable solution and
# the derivative is not too close to zero
while( steps<=iterations and abs(f(xn))>tolerance and abs(fprime(xn))>tolerance ):
xn -= f(xn)/fprime(xn)
myList += [(xn,0)]
steps += 1
#output the results
if(abs(fprime(xn))<= tolerance ):
print ("error, derivative value close to zero", xn.n())
else:
print ("approximate solution = ", xn.n())
#make the plot either way - so it can be analyzed
# - plot the list of points
p = list_plot(myList, color = "red", size = 50)
# - to put the text, we need to find how big is the graph
y1 = min([f(myList[i][0]) for i in [0,1,..,len(myList)-1]])
y2 = max([f(myList[i][0]) for i in [0,1,..,len(myList)-1]])
posT = y1-0.10*(y2-y1)
t = text(1, (myList[0][0], posT), fontsize=20, color = "black")
for i in range(1, len(myList)):
t = t+ text(i+1, (myList[i][0],posT), fontsize=20,
color = "black")
X1 = myList[i-1][0]
Y1 = f(X1)
X2 = myList[i][0]
Y2 = f(X2)
#tangent lin"
t = t + line([(X1, Y1),(X2,0)], color = "green")
#vertical lin"
t = t + line([(X2, 0),(X2,Y2)], color = "green")
# - plot the function
a = min([myList[i][0] for i in [0,1,..,len(myList)-1]])
b = max([myList[i][0] for i in [0,1,..,len(myList)-1]])
q = plot(f(x), x, a-5, b+5)
(p+q+t).show()
In [18]:
# Newtonova Raphsonova metóda na riešenie rovníc - iteračné kroky algoritmu
# Vstupy
f(x) = x^2-7
x1=30
iterations=1
tolerance=0.00001
#we start with x1
xn = x1
#this will collect the list of all successive approximations
myList = [(xn,0)]
#we need the derivative of f with respect to x
fprime(x) = derivative(f(x),x)
steps = 1
#keep computing successive approximations as long as
# we did not reach max_steps
# we do not have a reasonable solution and
# the derivative is not too close to zero
while( steps<=iterations and abs(f(xn))>tolerance and abs(fprime(xn))>tolerance ):
xn -= f(xn)/fprime(xn)
myList += [(xn,0)]
steps += 1
#output the results
if(abs(fprime(xn))<= tolerance ):
print ("error, derivative value close to zero", xn.n())
else:
print ("approximate solution = ", xn.n())
#make the plot either way - so it can be analyzed
# - plot the list of points
p = list_plot(myList, color = "red", size = 50)
# - to put the text, we need to find how big is the graph
y1 = min([f(myList[i][0]) for i in [0,1,..,len(myList)-1]])
y2 = max([f(myList[i][0]) for i in [0,1,..,len(myList)-1]])
posT = y1-0.10*(y2-y1)
t = text(1, (myList[0][0], posT), fontsize=20, color = "black")
for i in range(1, len(myList)):
t = t+ text(i+1, (myList[i][0],posT), fontsize=20,
color = "black")
X1 = myList[i-1][0]
Y1 = f(X1)
X2 = myList[i][0]
Y2 = f(X2)
#tangent lin"
t = t + line([(X1, Y1),(X2,0)], color = "green")
#vertical lin"
t = t + line([(X2, 0),(X2,Y2)], color = "green")
# - plot the function
a = min([myList[i][0] for i in [0,1,..,len(myList)-1]])
b = max([myList[i][0] for i in [0,1,..,len(myList)-1]])
q = plot(f(x), x, a-5, b+5)
(p+q+t).show()
