🧮 Calculus v jazyku Python

Knižnica SymPy

Balíček/knižnica SymPy vo všeobecnosti slúži na vykonávanie symbolických operácii v jazyku Python.

import sympy as sp # Import knižnice SymPy s premenovaním na skratkové slovo - sp

U,V = sp.symbols("U V") # Definovanie premenných U,V ako symbolické premenné z knižnice SymPy

Alternativne môžeme symbolické premenné definovať aj samostatne nasledovnou syntaxou:

U = sp.Symbol("U")
V = sp.Symbol("V")

Vo všeobecnosti je tak definovanie symbolických premenných možné opísať nasledovne:

NÁZOV_PREMENNEJ_1, NÁZOV_PREMENNEJ_2... = sp.symbols("NÁZOV_PREMENNEJ_1 NÁZOV_PREMENNEJ_2...") (Názvy premenných vo funkcii symbols v úvodzovkách sú oddelené medzerami)

Ak chceme jednotlivé symbolické premenné definovať samostatne je syntax nasledovná:

NÁZOV_PREMENNEJ = sp.Symbol("NÁZOV_PREMENNEJ")

Ďalej s premennými vo vzťahoch pracujeme ako štandardne. Avšak pozor na použitie ostatných matematických funckií. Tie musia byť taktiež symbolické. Ukážeme si to na priklade definície zobrazovacích rovníc z cvičenia 3.

from Degrees import *

P = {"U":np.radians(dms2deg(35,33,0)),"V":np.radians(dms2deg(15,8,0))} # Numerická premenná
R = 6371000 # Numerická premenná

x = ((R/2)*((sp.pi/2)-U))*(V+sp.sin(V)) # Definovanie zobrazovacej rovnice pre x
y = -((R/2)*((sp.pi/2)-U))*(1+sp.cos(V)) # Definovanie zobrazovacej rovnice pre y 

Všimnime si, že zo symbolickými premennými U a V pracujeme analogicky ako v štandardnom prípade kedy ide o numerické premenné.

V prípade použitia iných matematických funckcí alebo konštánt v tomto zápise je však taktiež nutné použiť symbolický typ. Konkrétne napr. pre výpočet kosínusu uhla by sme využili funkciu z knižnice NumPy, tento postup by v tomto prípade nefungoval. Funckie v symbolických výrazoch je tak potrebné využiť z knižnice SymPy, podobne ako to vydíte v rovniciach vyššie v prípade sínusu a konštanty pí.

---

Parciálne derivácie pomocou knižnice SymPy

Za predpokladu, že premenné máme definované symbolicky môžeme pristúpiť k derivácii výrazu. Na deriváciu funkcie podľa požadovanej premennej využijeme funkciu sp.diff(). Jej použitie je nasledovné:

sp.diff(DERIVOVANÁ_FUNKCIA, PREMENNÁ_PODĽA_KTOREJ_DERIVUJEME)

Prakticky si to môžeme ukázať na príklade z cvičenia, derivujem funkciu(výraz) x podľa premennej U:

dx_U = sp.diff(x,U)

Takto zderivovaná fukcia/výraz je však dzerivovaný len symbolický. Za predpokladu, že výsledkom derivácie nie je len konštanta, vo výsledku stále figurujú premenné. Pre ďalšiu prácu je tak potrebné určiť hodnotu derivácie dosadením numerických hodnôt. To vykonáme substitúciou(dosadením), pomocou metódy .subs(). Tento postup sa mierne líši od použitia bežných funkcií. Vzľadom na objektovo orientovaný prístup balíčka SymPy je potrebné substitúciu vykonať nasledovne:

KAM_CHCEM_DOSADIŤ.subs(ZA_ČO_DOSADZUJEM, HODNOTA_KTORÚ_DOSADZUJEM)

V prípade (je to aj náš prípad), že chcem dosadiť viacero hodnôt, je potrebné použiť zoznam(list) n-tíc(tuples), syntax je nasledovná:

KAM_CHCEM_DOSADIŤ.subs([(ZA_ČO_DOSADZUJEM_1, HODNOTA_KTORÚ_DOSADZUJEM),(ZA_ČO_DOSADZUJEM_2, HODNOTA_KTORÚ_DOSADZUJEM)...])

V našom prípade to urobíme nasledovne:

dx_U = dx_U.subs([(U,P["U"]),(V,P["V"])])

Takýmto spôsobom získame hodnotu parciálnej derivácie. Má to však jeden problém. Takýto výsledok má definovaný špecifický objektový dátový typ z balíčka SymPy. Pre ďalšie výpočny je tak nie veľmi dobre použiteľný. Opravíme to použitím funckie float(), ktorá z tohot dátového typu urobí štandardný dátový typ čísla z plávjúcou desatinnou čiarkou. Syntax je nasledovná: float(PREMENNÁ)

dx_U = float(dx_U)

Celý postup môžeme zapísať aj pomocou rekurentných definícií nasledovne:

dx_U = float(sp.diff(x,U).subs([(U,P["U"]),(V,P["V"])]))

---

Určité integrály pomocou knižnice SymPy

Podobne ako sme realizovali deriváciu pomocou knižnice SymPy tak podobným spôsobom vieme aj integravať, využijeme na to funckiu sp.integrate(). Syntax jej použitia je nasledovná:

sp.integrate(INTEGROVANÁ_FUNKCIA, (INTEGRAČNÁ_PREMENNÁ, DOLNÁ_HRANICA, HORNÁ_HRANICA))

Nezabudnite si všimnúť, že v prípade výpočtu určitého integrálu je druhý argument funkcie sp.integrate() definovaný ako n-tica(tuple) s potrebnými hodnotami, konkrétne

(INTEGRAČNÁ_PREMENNÁ, DOLNÁ_HRANICA, HORNÁ_HRANICA)

# Definovanie nutných vstupov pre ukážku
GRS80 = {"a":6378137, "e":0.08181919112}
Havana = {"phi":np.radians(dms2deg(23,8,0)), "lambda":np.radians(dms2deg(-82,23,0))}
Bratislava = {"phi":np.radians(dms2deg(48,9,0)), "lambda":np.radians(dms2deg(17,7,0))}
A_HB = np.radians(dms2deg(72,36,41.853))

Pre optimalizáciu som výpočet zakomentoval, žerie to dosť výpočtového výkonu.

phi = sp.Symbol("phi") # Symbolické definovanie premennej

#f = (GRS80["a"]*(1-(GRS80["e"]**2)))/(sp.sqrt((1-(GRS80["e"]**2)*sp.sin(phi)**2))**3) # Definovanie integrovanej funkcie

#l = (1/(np.cos(A_HB)))*(float(sp.integrate(f2,(phi,Havana["phi"],Bratislava["phi"]))))

Podobne, ako pri derivácii je vhodné použiť aj funkciu float(), avšak v tomto prípade to nie je vyložene nutné.

---

Určité integrály pomocou knižnice SciPy

Pre rýchlejší výpočet určitého integrálu je vhodné použiť funkciu z vedeckej knižnice SciPy. Syntax výpočtu je však mierne odlišná v porovnaní z výpočtami pomocou NumPy. Z tohto balíčka budeme potrebovať len jednu konkrétnu funkciu quad(). Preto je vhodné tento balíček importovať nasledovne:

from scipy.integrate import quad

K integrácii pristupíme tak, že najprv definujeme integrovanú matematickú funkciu, ako užívateľskú funkciu v jazyku python. Podrobnostiam tejto definície sa budeme venovať v nasledujúcej kapitole. Pre tento prípad však stačí využiť nasledujúcu syntax:

def NÁZOV_FUNKCIE(PREMENNÉ_ODDELENÉ_ČIARKOU):
    return ČO_MÁ_FUNKCIA_VRÁTIŤ

Kód pre náš prípad môžeme zapísať nasledovne:

def f(phi):
    return (GRS80["a"]*(1-(GRS80["e"]**2)))/(np.sqrt((1-(GRS80["e"]**2)*np.sin(phi)**2))**3)

Následne vieme na takto definovanú funckiu zavolať funkciu z balíčka SciPy quad(). Syntax je nasledovná:

quad(FUNKCIA,DOLNÁ_HRANICA,HORNÁ_HRANICA)

Výsledkom použitia tejto funkcie sú dve hodnoty. Pre ďalšie výpočty je tak nutné uvážiť správny index vypočítanej hodnoty. V našom prípade sú indexy nasledovné:

• O = Numerická hodnota určitého integrálu s dosadenými hranicami
• 1 = Presnosť numerického výpočtu

quad(FUNKCIA,DOLNÁ_HRANICA,HORNÁ_HRANICA)[INDEX]

Pre naše potreby použijeme: quad(FUNKCIA,DOLNÁ_HRANICA,HORNÁ_HRANICA)[0]

Syntaktický zápis tak vyzerá nasledovne:

l = (1/(np.cos(A_HB)))*(quad(f,Havana["phi"],Bratislava["phi"])[0])

---

Pokročilejšie funkcie balíčka SciPy

V rámci minulého tutoriálu sme sa už venovali jednoduchej numerickej integrácii pomociu funkcie quad() z balíčka SciPy. Základné predpoklady pre prácu s matematickými funckcionalitami knižnice SciPy ostávajú nezmenené. Predpokladom je, že funkcia ktorú integrujeme, derivujeme… je definovaná ako užívateľská funkcia v prostredí Python. Na pripomenutie tu je všeobecná syntax:

def NÁZOV_FUNKCIE(PREMENNÉ_ODDELENÉ_ČIARKOU):
    return ČO_MÁ_FUNKCIA_VRÁTIŤ

V tomto tutoriáli si priblížíme dve nové funkcie, a tými budú funkcie nquad() a newton(). Pre použitie je potrebné obe funkcie importovat:

from scipy.integrate import nquad
from scipy.optimize import newton

Dvojné určité integrály

Dvojný (určitý) integrál vypočítame použitím funckie nquad(). Jej použitie spočíva definovaní integrovanej funkcie. Keďže ide o dvojný (resp. n-tý) integrál, je potrebné aby funkcia ktorú definujeme bola závislá od dvoch (resp. n) premenných podla ktorých budeme integrovať. Zapísať to môžeme nasledovne:

def INTEGROVANÁ_FUNCKCIA(INTEGRAČNÉ_PREMENNÉ_ODDELENÉ_ČIARKOU):
    return ZÁPIS_FUNKCIE

Je nutné podotknúť, že itegračná funkcia musí byť závislá len od definovaných premenných, keďže k ostatným premenným sa bude pristupovať ako ku konštantám. Použitie môžeme inlustorvať na príklade s cvičení (výpočet plochy elipsoidického lichobežníka, ohraničeného šírkou a dĺžkou):

#Vyhľadom na obsenciu niektorých premenných je kód zakomentovaný

# def integrand(phi,lam):
#    return meridian(ell,phi)*lateral(ell,phi)*np.cos(phi)

Následne na zadefinovanú intregračnú funkciu aplikujeme funkciu nquad(). Jej všeobecná syntax je nasledovná:

nquad(INTEGROVANÁ_FUNKCIA,[[DOLNÁ_HRANICA_1, HORNÁ_HRANICA_1],[DOLNÁ_HRANICA_2, HORNÁ_HRANICA_2]])

Pri použití na konkrétnej situácii je syntax nasledovná:

#Vyhľadom na obsenciu niektorých premenných je kód zakomentovaný

#area = nquad(integrand,[[phi_min,phi_max],[lambda_min,lambda_max]]) 

Newtonova metóda

Funkcia newton() implementuje Newtonovu metódu a umožňuje riešiť nelineárne rovnice. Ak je derivácia ťažko vypočítateľná, môže použiť numerické diferencovanie. Naše riešenie spočíva len vo využití numerického diferencovania, aj napriek tomu, že jedným vstupom funkcie newton() môže byť aj prvá analyticky určená derivácia. Funkciu, ktorá je vstupom pre Newtonovu metódu je potrebné uviesť v tvare \(f(x)=0\). Syntax je tak nasledovná: `newton(FUNKCIA,POČIATOČNÝ_ODHAD_KOREŇA) Použitie je tiež dobre ukázateľné aj na príklade z cvičení:

#Vyhľadom na obsenciu niektorých premenných je kód zakomentovaný

# def fx_U(U):
#     return n*(rho_0*((np.tan((U_0/2)+np.radians(45)))/(np.tan((U/2)+np.radians(45))))**n)-(R*np.cos(U))

# U_1 = newton(fx_U,U_J)