# Numerické integrovanie a numerické riešenie diferenciálnych rovníc
# Programové kódy: SageMath v.9.x(Python 3.0)
# Prerekvizity: Integrálny počet, diferenciálne rovnice, 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 

# Výpočet neurčitého integrálu, výpočet primitívnej funkcie 
# Použitie knižnice "indefinite_integral"
from sage.symbolic.integration.integral import indefinite_integral

# vstup : funkcia f
# výstup: funkcia F

f(x)=x^3-3*x^2+2*x
r=LatexExpr("\\int(")+latex(f(x))+LatexExpr(")dx=")
F=indefinite_integral(f(x),x)
show(r+latex(F)+" + C")

$$\newcommand{\Bold}[1]{\mathbf{#1}}\int( x^{3} - 3 \, x^{2} + 2 \, x )dx= \frac{1}{4} \, x^{4} - x^{3} + x^{2} + C$$

# Výpočet neurčitého integrálu - vstup ľubovoľná funkcia cez @interact
# Tvorba tabuliek s integrálmi
# Použitie knižnice "indefinite_integral"
from sage.symbolic.integration.integral import indefinite_integral

# vstup : ľubovoľná spojitá funkcia f
# výstup: funkcia F

@interact
def integral(
    f = input_box(default = x^3-3*x^2+2*x),auto_update=False):
    f(x)=f
    r=LatexExpr("\\int(")+latex(f(x))+LatexExpr(")dx=")
    F=indefinite_integral(f(x),x)
    show(r+latex(F)+" + C")

# Výpočet určitého integrálu 
# Použitie knižnice "definite_integral"
# Primitívna funkcia, Newton - Leibnitzov vzorec, F(b)-F(a)

from sage.symbolic.integration.integral import definite_integral

# vstup : funkcia p a interval [a,b]
p=sin(x)*sin(x)
a=0
b=pi

r=LatexExpr("\\int_0^{\\pi}sin^2(x)dx=")
# Výpočet
q=definite_integral(p,x,a,b)
# Primitívna funkcia
F=diff(p,x)
show(r+"["+latex(F)+"]_0^{\\pi}="+latex(q))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\int_0^{\pi}sin^2(x)dx= [ 2 \, \cos\left(x\right) \sin\left(x\right) ]_0^{\pi}= \frac{1}{2} \, \pi$$

# Výpočet určitého integrálu funkcie f(x) 
# Niekedy je potrebné vykonať výpočet za určitých predpokladov
# Predpoklady môžeme nastaviť pomocou príkazu assume()
# Testovať platných predpokladov pomocou assumptions() 
# V riešení  erf - je tzn. error function

from sage.symbolic.integration.integral import definite_integral

var('x','a','b','c','y')
f(x)= e^(-pi*(x-2*y/c)^2/(a*a))
assume(a>0)
assume(c>0)
g(x) = definite_integral(f(x),x,a,b)
g.show()
assumptions()

$$\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{1}{2} \, a \operatorname{erf}\left(-\frac{\sqrt{\pi} {\left(a c - 2 \, y\right)}}{a c}\right) - \frac{1}{2} \, a \operatorname{erf}\left(-\frac{\sqrt{\pi} {\left(b c - 2 \, y\right)}}{a c}\right)$$

[a > 0, c > 0]

# Výpočet určitého integrálu, graf funkcie f a riešenie 
# Použitie knižnice "definite_integral"
from sage.symbolic.integration.integral import definite_integral

# vstup : funkcia p a interval [a,b]
f=sin(x)*sin(x)
a=0
b=pi

r=LatexExpr("\\int_0^{\\pi}sin^2(x)dx=")
# Výpočet
q=definite_integral(f,x,a,b)
show(r+latex(q))
s=plot(f,x,a,b,gridlines='minor', title="$ f(x)="+str(f)+" $", \
       fillcolor = "yellow", fill = true)
s.show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\int_0^{\pi}sin^2(x)dx= \frac{1}{2} \, \pi$$

# Určitý integrál - interaktívny výpočet s grafom 

# Použitie knižnice "definite_integral"
from sage.symbolic.integration.integral import definite_integral

# vstupy : funkcia f a interval [a,b].
# outputs : hodnota určitého integrálu f na  [a,b], oblasť pod (nad) f
# author: P. Szabó, 12 March 2018

@interact
def integral(
    f = input_box(default = x^3-3*x^2+2*x) ,
    a = input_box(default =  -0.5) ,
    b = input_box(default = 2.5),
    gr = selector(values = ["Yes", "No"],
                 label = "Graph", default = "No" ), auto_update=False):
    show("Integral computation")
    f(x)=f
    r=LatexExpr("\\int_a^{b}(")+latex(f(x))+LatexExpr(")dx=")
    q=definite_integral(f(x),x,a,b)
    show(r+latex(q))
    if (gr == "Yes"):
        r=plot(f(x),x,a,b,gridlines='minor',title="$f(x)="+str(f(x))+"$",\
               fillcolor = "yellow", fill = true)
        r.show()

# Riemannov integrál, delenie intervalu, integrálne súčty, limita postupnosti
# integrálnych súčtov, pozri [1].
# Vstupy: funkcia f, interval [a,b], n - počet delení intervalu
# Výstup: graf funkcie s delením intervalu a integrálnymi súčtami 

f(x)=x^3-3*x^2+2*x
a=0
b=2
n=10

# kreslenie funkcie
p = plot(f, x, a,b, color = "red", thickness = 4, alpha = 0.5, \
         fillcolor = "blue", fill = true)
dx = (b-a)/n
# kreslenie obdlžníkov
for i in [1,2,..,n]:
#for the shading ...
    xbegin = a+(i-1)*dx
    xend   = a+i*dx
    p += polygon( [(xbegin,0),  (xend,0),               \
                       (xend, f(xend)), (xbegin,f(xend))],  \
                       color="yellow",alpha = 0.5, aspect_ratio=
                       "automatic")
    #these are optional for a better graph:
    p += line([(xend,0), (xend, f(xend))], color = "black")
    p += line([(xend, f(xend)), (xbegin,f(xend))], color = "black")
    p += line([(xbegin, f(xend)),(xbegin,0) ], color = "black")
p.show(title = "Riemannove súčty s pravými koncovými bodmi", ticks = dx)

# Riemannov integrál, delenie intervalu, integrálne súčty, limita postupnosti
# integrálnych súčtov, pozri [1].
# Vstupy: funkcia f, interval [a,b], n - počet delení intervalu
# Vstupy interaktívne modifikovateľne
# Výstup: graf funkcie s delením intervalu a integrálnymi súčtami 

@interact
def RiemannSum(
    f = input_box(default = x^3-3*x^2+2*x) ,
    a = input_box(default =  0) ,
    b = input_box(default = 2) ,
    n = slider(vmin = 1, vmax = 20, default = 8, step_size=1)  ):
    f(x)=f
    #plot the function
    p = plot(f, x, a,b, color = "red", thickness = 4, \
             alpha = 0.5, fillcolor = "blue", fill = true)
    dx = (b-a)/n
    #plot each trapezoid
    for i in [1,2,..,n]:
    #for the shading ...
        xbegin = a+(i-1)*dx
        xend   = a+i*dx
        p += polygon( [(xbegin,0),  (xend,0),               \
                       (xend, f(xend)), (xbegin,f(xend))],  \
                       color="yellow",alpha = 0.5, aspect_ratio=
                       "automatic")
       #these are optional for a better graph:
        p += line([(xend,0), (xend, f(xend))], color = "black")
        p += line([(xend, f(xend)), (xbegin,f(xend))], color = "black")
        p += line([(xbegin, f(xend)),(xbegin,0) ], color = "black")
    p.show(title = "Riemann Sum with Right Endpoints", ticks = dx)

# Aplikácia určitého integrálu - výpočet objemu rotačného telesa
# Vstupy: funkcie f,g interval [a,b]
# Výstup: graf rotačného telesa medzi funkciami f,g na intervale
# [a,b], objem rotačného telesa
from sage.symbolic.integration.integral import definite_integral

f(x)=sqrt(2-x^2)
g(x)=x
a=0
b=1
# plot 2 funkcie
P1 =plot((f,g), (x,a,b), fill = True)
G2 = revolution_plot3d(f,(x,a,b), show_curve =True, \
            color='blue', opacity =0.5 ,parallel_axis ='x')
G3 = revolution_plot3d(g,(x,a,b), show_curve =True, \
            color = 'green', opacity =0.5 ,parallel_axis ='x')
G4 = G2+G3
# výpočet integrálu
Objem = pi*definite_integral((f(x)^2 - g(x)^2),x,a,b)
G4.show()
q=LatexExpr("Objem =")
show(q,latex(Objem))

$$\newcommand{\Bold}[1]{\mathbf{#1}}Objem = \frac{4}{3} \, \pi$$

# Aplikácia určitého integrálu - výpočet objemu rotačného telesa
# Vstupy: funkcie f,g interval [a,b] - interaktívne modifikovateľne
# Výstup: graf rotačného telesa medzi funkciami f,g na intervale
# [a,b], objem rotačného telesa
from sage.symbolic.integration.integral import definite_integral

@interact
def RotTeleso(
    f = input_box(default = x),
    g = input_box(default = x^2),
    a = input_box(default = 0),
    b = input_box(default = 1)):
    f(x)=f
    g(x)=g
    # plot 2 funkcie
    P1 =plot((f,g), (x,a,b), fill = True)
    G2 = revolution_plot3d(f,(x,a,b),show_curve =True, \
                color='blue', opacity =0.5 ,parallel_axis ='x')
    G3 = revolution_plot3d(g,(x,a,b),show_curve =True, \
                color = 'green', opacity =0.5 ,parallel_axis ='x')
    G4 = G2+G3
    # výpočet integrálu
    Objem = pi*definite_integral((f(x)^2 - g(x)^2),x,a,b)
    G4.show()
    print ('Objem =',Objem)

# Numerický výpočet určitého integrálu
# Lichobežníková metóda, delenie intervalu, pozri [1] 
# Odhad - súčet plošných obsahov lichobežníkov
# Používanie knižnice "numpy" - numerické výpočty pre Python
# Vstupy: funkcia, f, interval [a,b], n - počet delení intervalu

import numpy
f(x)=x**3*sin(x)
a=0
b=3.14
n=8

#compute the size of each subinterval
delta_x = (b-a)/n
#compute the intermediate points
xi_values = [a,a+delta_x,..,b]
#compute the corresponding y values
fxi_values = [f(i) for i in xi_values]
#compute the Trapezoidal Rule value
print ("approximation = ",  numpy.trapz(fxi_values ,
    dx=delta_x).n(), \
          "\t exact value = ",  integrate(f(x), x, a, b).n())
p = plot(f, x, a,b, color = "red", thickness = 2)    
#plot each trapezoid
p += line( [(a,0),  (a,f(a))] , color = "black")
for i in [1,2,..,n]:
    p += line( [(a+(i-1)*delta_x,f(a+(i-1)*delta_x)), \
        (a+i*delta_x,f(a+i*delta_x))] , color = "black")
    p += line( [(a+i*delta_x,0),  (a+i*delta_x,f(a+i*delta_x))] , \
        color = "black")
        #for the shading ...
    p += polygon( [(a+(i-1)*delta_x,0),              \
        (a+(i-1)*delta_x,f(a+(i-1)*delta_x)),   \
        (a    +i*delta_x, f(a+i*delta_x)),      \
        (a+    i*delta_x,       0)],      \
        color="yellow",             \
        aspect_ratio="automatic")
p.show()

approximation =  11.7606987535509 	 exact value =  12.1566814742107

# Numerický výpočet určitého integrálu
# Lichobežníková metóda, delenie intervalu, pozri [1]
# Odhad - súčet plošných obsahov lichobežníkov
# Používanie knižnice "numpy" - numerické výpočty pre Python
# Vstupy: funkcia, f, interval [a,b], n - počet delení intervalu
# Aplikácia @interact  - modifikácia vstupov
import numpy
@interact
def TrapezoidByNumPy(
    f = input_box(default = e^x),
    a = input_box(default = 0),
    b = input_box(default = 3),
    n = slider(vmin=2, vmax=50, default = 4, step_size = 1)):
    #needed to avoid the warning message
    f(x)=f
    #compute the size of each subinterval
    delta_x = (b-a)/n
    #compute the intermediate points
    xi_values = [a,a+delta_x,..,b]
    #compute the corresponding y values
    fxi_values = [f(i) for i in xi_values]
    #compute the Trapezoidal Rule value
    print ("approximation = ",  numpy.trapz(fxi_values ,
    dx=delta_x).n(), \
          "\t exact value = ",  integrate(f(x), x, a, b).n())
    p = plot(f, x, a,b, color = "red", thickness = 4)    
    #plot each trapezoid
    p += line( [(a,0),  (a,f(a))] , color = "black")
    for i in [1,2,..,n]:
        p += line( [(a+(i-1)*delta_x,f(a+(i-1)*delta_x)), \
                    (a+i*delta_x,f(a+i*delta_x))] , color = "black")
        p += line( [(a+i*delta_x,0),  (a+i*delta_x,f(a+i*delta_x))] , \
                   color = "black")
        #for the shading ...
        p += polygon( [(a+(i-1)*delta_x,0),              \
                       (a+(i-1)*delta_x,f(a+(i-1)*delta_x)),   \
                       (a    +i*delta_x, f(a+i*delta_x)),      \
                       (a+    i*delta_x,       0)],      \
                       color="yellow",             \
                       aspect_ratio="automatic")
    p.show()

# Symbolický výpočet určitého integrálu na intervale <0,oo)
u = var('u'); 
f = x * cos(u) / (u^2 + x^2);
assume(x>0); 
h = f.integrate(u, 0, infinity);
show(h)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \pi e^{\left(-x\right)}$$

# Numerický výpočet integrálu sin(x)/x, od 0 do 1
integral_numerical(sin(x)/x, 0, 1)

(0.946083070367183, 1.0503632079297087e-14)

# Integrál funkcie exp(-x**2), na intervale <0,oo)
# Symbolické a číselné vyjadrenie výsledku 
g = integrate(exp(-x**2), x, 0, infinity);
show(g) 
g.n()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \sqrt{\pi}$$

0.886226925452758

# Riešenie obyčajnej diferenciálnej rovnice
# vstupy: x'(t)+ x(t) = 2; pozri [1]
# výstup: funkcia x(t)
var('t')
x = function('x')(t)
de = diff(x,t) + x - 2
p=desolve(de==0, x)
q=LatexExpr("x(t)=")
show(q,p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}x(t)= {\left(C + 2 \, e^{t}\right)} e^{\left(-t\right)}$$

# Riešenie obyčajnej diferenciálnej rovnice
# vstupy: x'+ x = 2; pozri [1]
# výstup: funkcia x(t)
# Aplikácia @interact  - modifikácia vstupov
var('t')
@interact
def Interact2( de = input_box(default = "diff(x,t) + x - 2" ,
            label = "ODE = ")):
    x = function('x')(t)
    desolve(de==0, x).show()

# Riešenie obyčajnej diferenciálnej rovnice
# vstupy: x'+ x = 2; x(0) = 10 (začiatočná podmienka)
# výstup: funkcia x(t)
var('t')
x = function('x')(t)
de = diff(x,t) + x - 2
p=desolve(de==0, x, ics = [0,10])
q=LatexExpr("x(t)=")
show(q,p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}x(t)= 2 \, {\left(e^{t} + 4\right)} e^{\left(-t\right)}$$

# Riešenie obyčajnej diferenciálnej rovnice
# vstupy: x'+ x = 2; x(0) = 10 (začiatočná podmienka)
# výstup: funkcia x(t)
# aplikácia @interact - interaktívne vstupy, pozri [1] 

var('t')
x = function('x')(t)
@interact
def ODEInteract(
    choice = selector(values = ['yes','no'],
                      label = 'initial condition?',
                      default = 'yes')):
    if (choice == 'yes'):
        @interact
        def Interact1( de = input_box(default = "diff(x,t) + x - 2",
                                      label = "ODE = "),
                       t0 = input_box(default = 0,
                                       label = "t0="),
                       xt0 = input_box(default = 10,
                                       label = "x(t0)=")):
             desolve(de==0, x, ics = [t0, xt0]).show()
    else:
        @interact
        def Interact2( de = input_box(default = "diff(x,t) + x - 2",
                label = "ODE = ")):
            desolve(de==0, x).show()

# Nasledujúce príklady sa nachádzajú v knihe 
# Computational Mathematics with SageMath,  pozri[3]
# Základné parametre a  kód systému riešenia dif. rovníc
# x = var('x')
# y = function('y')(x)
# desolve(equation, variable, ics = ..., ivar = ..., 
# .....show_method = ..., contrib_ode = ...)
# Lineárne diferenciálne rovnice typu: y' + P(x)y = Q(x) 
# Riešte dif. rovnicu y' + 3y = e^x
x = var('x') 
y = function('y')(x)
p = desolve(diff(y,x) + 3*y == exp(x), y, show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{4} \, {\left(4 \, C + e^{\left(4 \, x\right)}\right)} e^{\left(-3 \, x\right)}, \verb|linear|\right]$$

# Separovateľné diferenciálne rovnice typu: P(x) = y'Q(y) 
# Riešte dif. rovnicu x=y'y
x = var('x') 
y = function('y')(x)
p = desolve(y*diff(y,x) == x, y, show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, y\left(x\right)^{2} = \frac{1}{2} \, x^{2} + C, \verb|separable|\right]$$

# Bernoulliho diferenciálne rovnice typu: y' + P(x)y = Q(x)y^alpha, 
# alpha<>0,alpha<>1. 
# Riešte dif. rovnicu y' − y = xy^4
x = var('x') 
y = function('y')(x)
p = desolve(diff(y,x)-y == x*y^4, y, show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{e^{x}}{{\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} + C\right)}^{\frac{1}{3}}}, \verb|bernoulli|\right]$$

# Homogénne diferenciálne rovnice typu: y'= P(x,y)/Q(x,y) 
# Riešte dif. rovnicu x^2y' = y^2 + xy + x^2.
x = var('x') 
y = function('y')(x)
p = desolve(x^2*diff(y,x) == y^2+x*y+x^2, y, show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[C x = e^{\arctan\left(\frac{y\left(x\right)}{x}\right)}, \verb|homogeneous|\right]$$

# Exaktné diferenciálne rovnice typu:  ∂f/∂x dx + ∂f/∂y dy,
# Riešte dif. rovnicu y' = (cos(y)−2x)/(y+x sin(y)),
# pričom f(x,y)= x^2 − xcos(y) + y^2/2.
x = var('x') 
y = function('y')(x)
p = desolve(diff(y,x)==(cos(y)-2*x)/(y+x*sin(y)), y,show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x^{2} - x \cos\left(y\left(x\right)\right) + \frac{1}{2} \, y\left(x\right)^{2} = C, \verb|exact|\right]$$

# Zložitejšie diferenciálne rovnice riešíme tak, že nastavíme
# parameter  contrib_ode=True
# y' = P(x)y^2 + Q(x)y + R(x), (dif. rovnice typu Ricatti)
# vo výstupe sa nachádzajú tzn. Besselove funkcie
# Riešte dif. rovnicu:  y' = xy^2 +(1/x)y − 1/x^2 
x = var('x') 
y = function('y')(x)
p=desolve(diff(y,x) == x*y^2+y/x-1/x^2, y,contrib_ode=True,\
          show_method=True)
print(p)
show(p)

[[y(x) == -1/2*((_C*(bessel_Y(4, 2*sqrt(-x)) - bessel_Y(2, 2*sqrt(-x))) + bessel_J(4, 2*sqrt(-x)) - bessel_J(2, 2*sqrt(-x)))*x + 3*(_C*bessel_Y(3, 2*sqrt(-x)) + bessel_J(3, 2*sqrt(-x)))*sqrt(-x))/((_C*bessel_Y(3, 2*sqrt(-x)) + bessel_J(3, 2*sqrt(-x)))*sqrt(-x)*x^2)], 'riccati']

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left[y\left(x\right) = -\frac{{\left(C {\left(Y_{4}(2 \, \sqrt{-x}) - Y_{2}(2 \, \sqrt{-x})\right)} + J_{4}(2 \, \sqrt{-x}) - J_{2}(2 \, \sqrt{-x})\right)} x + 3 \, {\left(C Y_{3}(2 \, \sqrt{-x}) + J_{3}(2 \, \sqrt{-x})\right)} \sqrt{-x}}{2 \, {\left(C Y_{3}(2 \, \sqrt{-x}) + J_{3}(2 \, \sqrt{-x})\right)} \sqrt{-x} x^{2}}\right], \verb|riccati|\right]$$

# Lineárne diferenciálne rovnice  
# Riešte dif. rovnicu y' + 2y = x^2 − 2x + 3:
x = var('x') 
y = function('y')(x)
DE = diff(y,x)+2*y == x**2-2*x+3
p = desolve(DE, y,show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{4} \, {\left({\left(2 \, x^{2} - 2 \, x + 1\right)} e^{\left(2 \, x\right)} - 2 \, {\left(2 \, x - 1\right)} e^{\left(2 \, x\right)} + 4 \, C + 6 \, e^{\left(2 \, x\right)}\right)} e^{\left(-2 \, x\right)}, \verb|linear|\right]$$

# Lineárne diferenciálne rovnice  
# Riešte dif. rovnicu y' + 2y = x^2 − 2x + 3:
# zjednodušenie riešenia
x = var('x') 
y = function('y')(x)
DE = diff(y,x)+2*y == x**2-2*x+3
p = desolve(DE, y).expand()
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, x^{2} + C e^{\left(-2 \, x\right)} - \frac{3}{2} \, x + \frac{9}{4}$$

# Riešte dif. rovnicu y'log(y) = y sin(x)
x = var('x')
y = function('y')(x)
de = diff(y,x)*log(y) == y*sin(x)
p = desolve(de, y, show_method=True)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{2} \, \log\left(y\left(x\right)\right)^{2} = C - \cos\left(x\right), \verb|separable|\right]$$

# Riešte dif. rovnicu y'log(y) = y sin(x)
# Riešenie úlohy sa nachádza v zloženej funkcii
# Preto ešte treba vyjadriť riešenie, funkciu y(x)
x = var('x')
y = function('y')(x)
ed = desolve(diff(y,x)*log(y) == y*sin(x), y)
p = solve(ed,y)
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[y\left(x\right) = e^{\left(-\sqrt{2 \, C - 2 \, \cos\left(x\right)}\right)}, y\left(x\right) = e^{\left(\sqrt{2 \, C - 2 \, \cos\left(x\right)}\right)}\right]$$

# Riešenie je definované pre C>1. 
# Nasledovný kód nakreslí riešenie pre C=2
# Najprv zistíme premenné (konštanty) v riešení
ed.variables()

(_C, x)

# Teraz dosadíme do konštanty c, c=5
c = ed.variables()[0]
p=solve(ed, y)[0].substitute(c == 5).rhs()
show(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}e^{\left(-\sqrt{-2 \, \cos\left(x\right) + 10}\right)}$$

# Nakreslenie riešenia ODE pre niektoré konštanty c.

p = plot(solve(ed, y)[0].substitute(c == 5).rhs(), x, -3, 3)
p +=plot(solve(ed, y)[0].substitute(c == 4.5).rhs(), x, -3, 3)
p +=plot(solve(ed, y)[0].substitute(c == 4).rhs(), x, -3, 3)
show(p)

# Eulerova metóda na numerické riešenie obyčajných diferenciálnych rovníc
# Vstupy: funkcia y'=f(x,y), interval [a,b], začiatočná hodnota [a,y(a)]
# n - počet bodov v tabuľke pre funkciu y(x), pozn: y(a)=y(0)
# xk= y(xk-1)+h*f(xk-1,yk-1); h=(b-a)/n
# výstup: tabuľka bodov funkcie y(x): [[x0,y0],[x1,y1],..,[xn,yn]]
# 29.11.2018, P. Szabó
var("x,y")
f(x,y)=cos(x)+sin(y)+x*y
n=15
a=0
b=2
y0=8 
#compute the step size
h = (b-a)/n
#create the table of approximations
table = [(a, y0)]
prevx = a
prevy = y0
for i in [1,2,..,n]:
    #compute the new values
    newx = prevx+h
    newy = prevy + h*f(prevx, prevy)
    #add them the to the table
    table += [(newx, newy)]
    #"move" to the new values
    prevx = newx.n()
    prevy = newy.n()
    #plot the obtained table
list_plot(table,size = 50).show()

# Eulerova metóda na numerické riešenie obyčajných diferenciálnych rovníc
# Vstupy: funkcia y'=f(x,y), interval [a,b], začiatočná podmienka [a,y(a)]
# n - počet bodov v tabuľke pre funkciu y(x), pozn: y(a)=y(0)
# výstup: tabuľka bodov funkcie y(x): [[x0,y0],[x1,y1],..,[xn,yn]]
# aplikácia @interact - interaktívne vstupy, pozri[1]
# 29.11.2018, P. Szabó
var("x,y")
@interact
def EulerDEInteract(
    f = input_box(default= cos(x)+sin(y)+x*y),
    n = slider(vmin=0, vmax=100, step_size=1,
               default=3,
               label="Select the order n: "),
    a = input_box(default= 0),
    b = input_box(default = 2),
    y0= input_box(default = 8)):
    #this it to avoid a warning message
    f(x,y)=f
    #compute the step size
    h = (b-a)/n
    #create the table of approximations
    table = [(a, y0)]
    prevx = a
    prevy = y0
    for i in [1,2,..,n]:
        #compute the new values
        newx = prevx+h
        newy = prevy + h*f(prevx, prevy)
        #add them the to the table
        table += [(newx, newy)]
        #"move" to the new values
        prevx = newx.n()
        prevy = newy.n()
    #plot the obtained table
    list_plot(table,size = 50).show()

# Štruktúrované riešenie diferenciálnych rovníc
# T - je objekt na riešenie obyčajných diferenciálnych rovníc (ODE)
# T.function - je funkcia na pravej strane
# T.algorithm - je algoritmus na riešenie ODE "rk8pd"
#                Runge Kutta Prince Domand method
#           štandardná metóda na numerické riešenie ODE (2022 dec.)
# T.ode_solve - riešenie, iteratívne nájdenie funkčných hodnôt funkcie
# T.interpolate_solution - interpolačná funkcia podľa nájdených bodov
T = ode_solver()
def f_1(t,y,params): return [y[1],params[0]*(1-y[0]^2)*y[1]-y[0]]
T.function = f_1
# def j_1(t,y,params):
#    return [[0, 1], [-2*params[0]*y[0]*y[1]-1, params[0]*(1-y[0]^2)],[0,0]]
# T.jacobian = j_1
T.algorithm = "rk8pd"
T.ode_solve(y_0=[1,0], t_span=[0,100], params=[10],num_points=1000)
f = T.interpolate_solution()
plot(f, 0, 100)

# Interpolácia - Lagrangeov polynóm
# Na vstupe máme hodnoty [xi,yi], i=0,1,,.,n-1
# Na výstupe polynóm f rádu n-1 pre ktorý f(xi)=yi, i=0,1,,.,n-1
# Y pohľadu programovania je to ukážka ako aplikovať zabudované funkcie
# Vstupné body je možné modifikovať  - @interact

def Lagrange_Basis(i, points):
    var("x")
    n = len(points)-1
    Li=1
    for j in [0,1,..,n] :
        if(i!=j):
            Li *= (x-points[j][0])/(points[i][0]-points[j][0])
    return Li

def Lagrange_Polynomial(points):
    var("x")
    n = len(points)-1
    p=0*x
    for i in [0,1,..,n]:
        p +=  points[i][1]*Lagrange_Basis(i,points)
    return p.full_simplify()

@interact
def LagrangeInterpolationInteract(
    points = input_box(default = [(-3,-15),(-1,-5),(0,1),(2,10),
    (3,15)]) ):
    # print "p(x)=", Lagrange_Polynomial(points)
    r2=LatexExpr(" f(x)=")+latex(Lagrange_Polynomial(points))
    show(r2)

# Interpolácia - Lagrangeov polynóm 
# (existujú aj iné interpolačné algoritmy)
# Na vstupe máme hodnoty [xi,yi], i=0,1,,.,n-1
# Na výstupe polynóm f rádu n-1 pre ktorý f(xi)=yi, i=0,1,,.,n-1

def Lagrange_Basis(i, points):
    var("x")
    n = len(points)-1
    Li=1
    for j in [0,1,..,n] :
        if(i!=j):
            Li *= (x-points[j][0])/(points[i][0]-points[j][0])
    return Li

def Lagrange_Polynomial(points):
    var("x")
    n = len(points)-1
    p=0*x
    for i in [0,1,..,n]:
        p +=  points[i][1]*Lagrange_Basis(i,points)
    return p.full_simplify()


points = [(-3,-15),(-1,-5),(0,1),(2,10), (3,15)]
# print "p(x)=", Lagrange_Polynomial(points)
r2=LatexExpr(" f(x)=")+latex(Lagrange_Polynomial(points))
show(r2)

$$\newcommand{\Bold}[1]{\mathbf{#1}} f(x)= \frac{1}{18} \, x^{4} - \frac{1}{18} \, x^{3} - \frac{11}{18} \, x^{2} + \frac{11}{2} \, x + 1$$

# Laplaceova transformácia
# vstup: funkcia f(t)
# výstup: Laplaceov obraz funkcie f(t)

f(t) = t*cos(t)
L(s) = laplace(f, t, s)
q=LatexExpr('\\mathcal{L}\\{t*cos(t)\\}=')
show(q,L(s))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathcal{L}\{t*cos(t)\}= \frac{2 \, s^{2}}{{\left(s^{2} + 1\right)}^{2}} - \frac{1}{s^{2} + 1}$$

# Laplaceova transformácia
# vstup: F(s) Laplaceov obraz funkcie 
# výstup: f(t) Laplaceov originál

L(s) = 5*e^(-3)/(s - 1)^2
f(t) = inverse_laplace( L(s), s, t )
q=LatexExpr('\\mathcal{L_{-1}}\\{5*e^{-3}/(s - 1)^2\\}=')
show(q,f(t))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathcal{L_{-1}}\{5*e^{-3}/(s - 1)^2\}= 5 \, t e^{\left(t - 3\right)}$$