#!/usr/bin/env python
# coding: utf-8

# <div style="font-size:35px"> Greške kod računanja vrijednosti funkcije    </div><br>

# 
# Neka je zadana funkcija $n$ promjenljivih $f(x_1,\ldots,x_n)$ definisana i diferencijabilna po svim argumenima u nekoj oblasti. Pretpostavimo da nisu poznate tačne vrijednosti promjenljivih $x_1,\ldots,x_n,$ nego umjesto njih znamo približne vrijednosti $x^{\star}_1,\ldots,x^{\star}_n,$ takođe smatrajmo da znamo i vrijednost greške približnog broja $x_i-x^{\star}_i$ i apsolutne greške $\triangle x^{\star}=\left|x_i-x^{\star}_i\right|,\:i=1,\ldots,n.$  
# Potrebno je izračunati vrijednost apsolutne greške funkcije $f$ u tački $x^{\star}=(x^{\star}_1,\ldots,x^{\star}_n),$ gdje je $x^{\star}_i,$  iz relacije   $x_i=x_i^{\star}\pm\triangle x_i^{\star},\,i=1,\ldots,n.$ 
# 
# Na osnovu Lagrangeove teoreme o srednjoj vrijednosti vrijedi
# $$
# \triangle f^{\star}=\left|f(x_1,\ldots,x_n)-f(x^{\star}_1,\ldots,x^{\star}_n)\right|
#      =\left|\sum_{i=1}^{n}{\frac{\partial f\left(\tilde{x}_1,\ldots,\tilde{x}_n \right)}{\partial x_i}}(x_i-x^{\star}_i)\right|,
# $$
# gdje je $\tilde{x}_i$ između $x_i$ i $x^{\star}_i.$ Dalje vrijedi
# $$
#  \left|\sum_{i=1}^{n}{\frac{\partial f\left(\tilde{x}_1,\ldots,\tilde{x}_n \right)}{\partial x_i}}(x_i-x^{\star}_i) \right|
#  \leqslant \sum_{i=1}^{n}{\left|\frac{\partial f\left(\tilde{x}_1,\ldots,\tilde{x}_n \right)}{\partial x_i}\right|}\left|x_i-x^{\star}_i\right|
#  =\sum_{i=1}^{n}{\left|\frac{\partial f\left(\tilde{x}_1,\ldots,\tilde{x}_n \right)}{\partial x_i}\right|}\triangle x^{\star}_i
# $$
# Pošto vrijednosti $\tilde{x}_i$ nisu poznate, umjesto $\frac{\partial f\left(\tilde{x}_1,\ldots,\tilde{x}_n \right)}{\partial x_i},$ koristićemo 
# $\frac{\partial f\left(x^{\star}_1,\ldots,x^{\star}_n \right)}{\partial x_i},$ pa dobijamo formulu za apsolutnu grešku vrijednosti funkcije
# $$
# \triangle f^{\star}\approx\sum_{i=1}^{n}{\left|\frac{\partial f\left(x^{\star}_1,\ldots,x^{\star}_n \right)}{\partial x_i}\right|}\triangle x^{\star}_i,(\star)
# $$
# a relativnu grešku približne vrijednosti funkcije je sada
# $$
#  \delta f^{\star}\approx\frac{\triangle f^{\star}}{\left| f^{\star}\right|},(\star\star)
# $$
# gdje je $f^{\star}=f(x^{\star}_1,\ldots,x^{\star}_n).$

# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA;padding: 5px 5px 5px 15px;" >    
#  <b>Napomena. </b> 
# Umesto relacija datih u $(\star)$ i $(\star \star)$, mogu se koristiti i sljedeće relacije   
# $$
# \triangle f^{\star}=\sum_{i=1}^{n}{\left|\frac{\partial f\left(x^{\star}_1,\ldots,x^{\star}_n \right)}{\partial x_i}\right|}\triangle x^{\star}_i,
# $$
# i 
# $$
# \delta f^{\star}=\frac{\triangle f^{\star}}{|f^{\star}|}
# $$
# Ovakve greške se nazivaju linearna ocjena apsolutne i relativne greške vrijednosti funkcije, vidjeti više u npr. Bahvalov $[Bah77]$.
# </div>

# <div style=" border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA;padding: 5px 5px 5px 15px;">    
# <b> Zadatak.</b> Izračunati apsolutnu i relativnu grešku pri računanju zapremine lopte $V=\frac{4}{3}r^3\pi,$ ako je poluprečnik lopte $r=10.3\pm 0.02$ cm, a $\pi^{\star}=3.14.$ 
# </div>
# <br>
# <b>Rješenje.</b> Vrijedi 
# \begin{align*}
# & r^{\star}=10.3,\, \triangle r^{\star}=0.02,\, \pi^{\star}=3.14,\, \triangle \pi^{\star}\approx 0.00159265\\
# &   \frac{\partial V(r^{\star},\pi^{\star})}{\partial r}=4(r^{\star})^2\pi^{\star}\\
# &   \frac{\partial V(r^{\star},\pi^{\star})}{\partial \pi}=\frac{4}{3}(r^{\star})^3\\
# &   \triangle V^{\star}\approx \left|\frac{\partial V(r^{\star},\pi^{\star})}{\partial r} \right|\triangle r^{\star}+
#               \left|\frac{\partial V(r^{\star},\pi^{\star})}{\partial \pi} \right|\triangle \pi^{\star}
# \end{align*}

# In[6]:


import numpy as np
import sympy #as sym

r, pi=sympy.symbols('r, pi')

deltar=0.02
deltapi=0.00159265

izvodr=sympy.diff(zapremina,r)
izvodpi=sympy.diff(zapremina,pi)

zapremina=(4/3)*r**3*pi
deltav=abs(izvodr)*deltar+abs(izvodpi)*deltapi

deltaV=deltav.subs([(r,10.3),(pi,3.14)]) 
zap=zapremina.subs([(r,10.3),(pi,3.14)]) 
print('Apsolutna greška je  %f,'%deltaV, ' približna vrijednost zapremine je %f.'%zap)


# Vrijedi $$V=\frac{4}{3}\left( r^3\right)^{\star}(\pi)^{\star}\pm \Delta V^{\star}=4574.88\pm 28.97.$$ 
# Kako je $\Delta V^{\star}\approx 28.97=0.2897\cdot 10^{2}\leqslant \frac{1}{2}\cdot 10^{3-2+1},$ sigurne su samo dvije cifre pa možemo zaokružiti $V^{\star}\approx 4600.$ I na kraju relativna greška je 
# $$\delta V=\frac{\Delta V^{\star}}{|V^{\star}|}\approx 0.006332\approx 0.6\%.$$

# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA;padding: 5px 5px 5px 15px;">   
# <b> Zadatak.</b> Po kosinusnoj teoremi, ako su za neki trougao poznate dužine stranica $a$ i $b,$ te ugao izmedju njih $\gamma,$ dužina treće stranice $c$ računa se koristeći formulu $$ c^2=a^2+b^2-2ab\cos \gamma.$$ Ako je za trougao $\Delta ABC$ poznato $a=10.2\pm 0.05$ cm, $b=7.5\pm 0.05$ cm i $\gamma=31.3^{0}\pm0.05^{0},$ odrediti apsolutnu i relativnu grešku pri računanju stranice $c.$  
# </div>
# 
# 
# <b> Rješenje.</b> Vrijednosti uglova iz stepena u radijane pretvaraju se po formuli $\psi_{rad}=\frac{\psi_{step}\cdot\pi}{180}.$ Dužinu stranice $c$ računamo po formuli
# 
# $$c=\sqrt{a^2+b^2-2ab\cos \gamma}.$$
# 
# 
# Vrijedi 
# \begin{align*}
# &a^{\star}= 10.2,\, \triangle a^{\star}=0.05,\, b^{\star}=7.5, \, \triangle b^{\star}=0.05,\, \gamma^{\star}=\frac{31.3\cdot \pi}{180}, \, \triangle\gamma^{\star}=\frac{0.05\cdot\pi}{180}  \\
# &\frac{\partial c(a^{\star},b^{\star}, \gamma^{\star})}{\partial a}
# =\frac{a^{\star}-b^{\star}\cos\gamma^{\star}}{\sqrt{(a^{\star})^2+(b^{\star})^2-2a^{\star}b^{\star}\cos\gamma^{\star}}}\\
# &\frac{\partial c(a^{\star},b^{\star}, \gamma^{\star})}{\partial b}
# =\frac{b^{\star}-a^{\star}\cos\gamma^{\star}}{\sqrt{(a^{\star})^2+(b^{\star})^2-2a^{\star}b^{\star}\cos\gamma^{\star}}}\\
# &\frac{\partial c(a^{\star},b^{\star}, \gamma^{\star})}{\partial \gamma}
# =\frac{a^{\star}b^{\star}\sin\gamma^{\star}}{\sqrt{(a^{\star})^2+(b^{\star})^2-2a^{\star}b^{\star}\cos\gamma^{\star}}}\\
# &\triangle c^{\star}=\left|\frac{\partial c(a^{\star},b^{\star}, \gamma^{\star})}{\partial a} \right|\triangle a^{\star}
#         +\left|\frac{\partial c(a^{\star},b^{\star}, \gamma^{\star})}{\partial b} \right|\triangle b^{\star}
#         +\left|\frac{\partial c(a^{\star},b^{\star}, \gamma^{\star})}{\partial \gamma} \right|\triangle \gamma^{\star}
# \end{align*}

# In[15]:


import numpy as np
from sympy import * #drugi način učitavanja biblioteke 

a,b,gama=symbols('a,b,gama')

stranicaC=sqrt(a**2+b**2-2*a*b*cos(gama))

deltaA=.05;deltaB=.05; deltaGama=(.05*np.pi)/180;

izvodA=diff(stranicaC,a)
izvodB=diff(stranicaC,b)
izvodGama=diff(stranicaC,gama)


deltac=abs(izvodA)*0.05+abs(izvodB)*deltaB+abs(izvodGama)*deltaGama

deltaC=deltac.subs([(a,10.2),(b,7.5),(gama,(31.3*np.pi/180))])
c=stranicaC.subs([(a,10.2),(b,7.5),(gama,(31.3*np.pi/180))])

print('Apsolutna greška je %f,'%deltaC, ' približna dužina stranice c je %f.'%c)


# Vrijedi $\Delta c^{\star}\approx0.052428=0.052428\cdot 10^{0}\leqslant 0.5\leqslant 0.5\cdot 10^{0}
#      =0.5\cdot 10^{0-1+1},$ dakle sigurna je samo prva cifra, pa je možemo zaokružiti $c^{\star}\approx 5.$ Relativna greška je $$\delta c^{\star}\approx \frac{\Delta c^{\star}}{|c^{\star}|}\approx 0.0096433\approx1\%.$$

# <div style="font-size:35px"> Inverzni (obratni) problem greške </div><br>
# 
# 
# 
# Inverzni problem greške je određivanje granice dopustivih vrijednosti grešaka argumenata za koje greška funkcije ne prelazi unaprijed zadanu vrijednost. U slučaju funkcija jedne realne promjenljive, ako je ta funkcija diferencijabilna vrijedi
# \begin{equation}
#  f(x)=f(x^{\star})+f'(\tilde{x})(x-x^{\star}), 
# \label{uvod24}
# \end{equation} 
# gdje je $\tilde{x}$ između $x$ i $x^{\star}.$ Sada za $f'(\tilde{x})\neq 0,$ vrijedi
# \begin{equation}
#  x-x^{\star}=\frac{f(x)- f(x^{\star})}{f'\left( \tilde{x}\right) }.
# \label{uvod25}
# \end{equation}  
# Kako ne znamo $\tilde{x},$ to je
# \begin{equation}
#  \triangle x^{\star}\approx\frac{\triangle f^{\star}}{\left|f'\left(x^{\star}\right) \right|},\:f'\left(x^{\star}\right) \neq 0.
# \label{uvod26}
# \end{equation}
# Ako je $f$ funkcija od $n$ argumenata, tj. $f(x_1,\ldots,x_n),$ onda se zadavanjem greške funkcije zadaje samo jedna veza između $n$ nepoznatih $\triangle x^{\star}_1,\ldots,\triangle x^{\star}_n.$  Ako je poznata apsolutna greška vrijednosti funkcije $\triangle f^{\star}$, dodatni uslovi koje apsolutne greške argumenata treba da zadovoljavaju obično se definišu na jedan od sljedećih načina
# 
# <br><b>(1) Princip jednakih uticaja (doprinosa, efekata)</b><br>
# Neka vrijedi
#    $$\left| \frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_1}\right|\triangle x^{\star}_1=
#             \cdots=\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_n}\right|\triangle x^{\star}_n,$$ te zbog $\triangle f^{\star}\approx\sum_{i=1}^{n}{\left|\frac{\partial f\left(x^{\star}_1,\ldots,x^{\star}_n \right)}{\partial x_i}\right|}\triangle x^{\star}_i,$  vrijedi 
#     $$\triangle f^{\star}\approx n\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_k} \right|\triangle x^{\star}_k,$$ i na kraju 
# $$
#        \triangle x^{\star}_k\approx\frac{\triangle f^{\star}}{n\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_k} \right|},\,k=1,\ldots,n.
# $$
#  <b> (2) Princip jednakih apsolutnih grešaka</b><br>
# Neka je sada 
#     $$\triangle x^{\star}_1=\cdots=\triangle x^{\star}_n,$$ ponovo zbog $\triangle f^{\star}\approx\sum_{i=1}^{n}{\left|\frac{\partial f\left(x^{\star}_1,\ldots,x^{\star}_n \right)}{\partial x_i}\right|}\triangle x^{\star}_i,$ imamo
# $$
#      \triangle f^{\star}\approx\triangle x^{\star}_k\sum_{j=1}^{n}{\left|\frac{\partial f(x^{\star}_1,  
#        \ldots,x^{\star}_n)}{\partial x_j} \right|},
# $$ te dobijamo
# $$
#        \triangle x^{\star}_k\approx\frac{\triangle f^{\star}}{\displaystyle\sum_{j=1}^{n}{\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_j} \right|}},\,k=1,\ldots,n.
# $$                
# <b> (3) Princip jednakih relativnih grešaka </b> <br>
# Neka je sada $$\delta x^{\star}_1=\cdots=\delta x^{\star}_n,$$ te ponovo zbog $\triangle f^{\star}\approx\sum_{i=1}^{n}{\left|\frac{\partial f\left(x^{\star}_1,\ldots,x^{\star}_n \right)}{\partial x_i}\right|}\triangle x^{\star}_i,$ imamo
# 
# \begin{align*}
# \triangle f^{\star}\approx&\sum_{j=1}^{n}{\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_j} \right|}
#                     \triangle x^{\star}_j\cdot \frac{\left|x^{\star}_j\right|}{\left|x^{\star}_j \right|}
#                   =\sum_{j=1}^{n}{\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_j} \right|}
#                     \left| x^{\star}_j\right|\cdot \frac{\triangle x^{\star}_j}{\left|x^{\star}_j \right|}\\
#                  =&\frac{\triangle x^{\star}_k}{\left|x^{\star}_k \right|}\cdot\sum_{j=1}^{n}{\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_j} \right|}
#                     \left| x^{\star}_j\right|,\,k=1,\ldots,n,   
# \end{align*}
# 
# pa je 
# \begin{equation*}
#  \triangle x^{\star}_k\approx\frac{\triangle f^{\star} \left|x^{\star}_k \right|}{\displaystyle\sum_{j=1}^{n}{\left|x^{\star}_j\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_j} \right|}},\,k=1,\ldots,n.
#  \label{uvod29}
# \end{equation*}
# 
# <div style=" border-left: 2px solid black; border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;    border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA; padding: 5px 5px 5px 15px">
# <b> Zadatak.</b>
# Potrebno je izradi pravilnu šestostranu piramidu, dužine stranice baze $a = 3.5$ cm i
# dužine visine bočne stranice $h = 7.2$ cm. Odrediti najveće dopuštene apsolutne greške promjenljivih $a$ i $h,$ tako da se površina izrad̄ene piramide razlikuje najviše za $\Delta P = 0.5$  od predvid̄ene. Zadatak uraditi koristeći sva tri principa.
# </div>
# 
# <b>Rješenje.</b> Površinu date šestostrane piramide računamo po formuli $$P=B+O=\frac{3\sqrt{3}a^2}{2}+3ah.$$
#  (1) Princip jednakih uticaja 
# 
# Vrijednost najveće dopuštene apsolutne greške argumenata računamo po formuli
# $$\triangle x^{\star}_k=\frac{\triangle f^{\star}}{n\left|\frac{\partial f{(x^{\star}_1,\ldots,x^{\star}_n)}}{\partial x_k} \right|},\,k=1,\ldots,n.$$ 
# U ovom slučaju to je
# 
# \begin{align*}
#   \triangle a^{\star}&\approx \frac{\triangle P^{\star}}{2\left| \frac{\partial P(a^{\star},h^{\star})}
#               {\partial a}\right|},\\
#   \triangle h^{\star}&\approx \frac{\triangle P^{\star}}{2\left| \frac{\partial P(a^{\star},h^{\star})}
#               {\partial h}\right|},\\
#   \frac{\partial P(a^{\star},h^{\star})}{\partial a}&=3\sqrt{3}a+3h,\\
#   \frac{\partial P(a^{\star},h^{\star})}{\partial h}&=3a.
# \end{align*}

# In[8]:


import sympy
import numpy as np

a,h=sympy.symbols('a,h')

deltaP=.5 

povrsina=1.5*np.sqrt(3)*a**2+3*a*h

izvoda=sympy.diff(povrsina,a)
izvodh=sympy.diff(povrsina,h)

deltaa=deltaP/(2*(abs( izvoda.subs([(a,3.5),(h,7.2)]))))
deltah=deltaP/(2*(abs( izvodh.subs([(a,3.5),(h,7.2)]))))

print('Najveća dopuštena vrijednost apsolutne greške stranice a je %1.7f,'%deltaa,
        '\nnajveća dopuštena vrijednost apsolutne greške visine h je %1.7f.'%deltah)


# (2) Princip jednakih apsolutnih grešaka  
# 
# Vrijednost najveće dopuštene apsolutne greške argumenata računamo po 
# $$\triangle x^{\star}_k
# =\frac{\triangle f^{\star}}{\sum_{j=1}^{n}\left|\frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n)}{\partial x_k} \right|   },\,k=1,\ldots,n,$$
# Pa je 
# $$\triangle a^{\star}=\triangle h^{\star}=\frac{\triangle P^{\star}}
#            {\left|\frac{\partial P(a^{\star}, h^{\star})}{\partial a} \right|
#                +\left|\frac{\partial P(a^{\star}, h^{\star})}{\partial h} \right|}.$$

# In[5]:


import sympy
import numpy as np

a,h =sympy.symbols('a,h')

deltaP=.5

povrsina=1.5*np.sqrt(3)*a**2+3*a*h

izvoda=sympy.diff(povrsina,a)
izvodh=sympy.diff(povrsina,h)

deltaah=deltaP/(abs(izvoda.subs([(a,3.5),(h,7.2)]))+abs(izvodh.subs([(a,3.5),(h,7.2)])))
print('Najveća dopuštena vrijednost apsolutne greške argumenata a i h je %1.7f.'%deltaah)


# (3) Princip jednakih relativnh grešaka 
# 
# Vrijednost najveće dopuštene apsolutne greške argumenata računamo po formuli 
# 
# $$\triangle x^{\star}_k=\frac{\triangle f^{\star}|x^{\star}_k|}
#  {\sum_{j=1}^{n}{\left|x^{\star}_j \frac{\partial f(x^{\star}_1,\ldots,x^{\star}_n )}{\partial x_j} \right|}},
#       \,k=1,\ldots,n, $$
# u ovom slučaju       
# \begin{align}
# \triangle a^{\star}&\approx\frac{\triangle P^{\star}|a^{\star}|}
#   {\left|a^{\star}\frac{\partial P(a^{\star},h^{\star})}{\partial a} \right|
#     +\left|h^{\star}\frac{\partial P(a^{\star},h^{\star})}{\partial h} \right|},\\
# \triangle h^{\star}&\approx\frac{\triangle P^{\star}|h^{\star}|}
#   {\left|a^{\star}\frac{\partial P(a^{\star},h^{\star})}{\partial a} \right|
#     +\left|h^{\star}\frac{\partial P(a^{\star},h^{\star})}{\partial h} \right|} .
# \end{align}

# In[17]:


import sympy
import numpy as np

a,h=sympy.symbols('a,h')
deltaP=.5;
povrsina=1.5*np.sqrt(3)*a**2+3*a*h

izvoda=sympy.diff(povrsina,a)
izvodh=sympy.diff(povrsina,h)

deltaa=(deltaP*a)/(abs(a*izvoda)+abs(h*izvodh))
deltah=(deltaP*h)/(abs(a*izvoda)+abs(h*izvodh))

print('Najveća dopuštena vrijednost apsolutne greške argumenta a je %1.7f,'%deltaa.subs([(a,3.5),(h,7.2)]),
      '\nnajveća dopuštena vrijednost apsolutne greške argumenta h je %1.7f'%deltah.subs([(a,3.5),(h,7.2)]))


# <div style="font-size:35px"> Literatura </div><br>
# 
# [Bah77]  N.S. Bahvalov, <i>Numerical Methods: Analysis, Algebra, Ordinary Differential Equations</i>, MIR, USSR, 1977. ISBN: 9780714712079.<br><br>
# [Sci15] R. Scitovski, <i>Numerička matematika,</i> Sveučilišta Josipa Jurja Strossmayera u Osijeku – Odjel za matematiku, R.Hrvatska, 2015. ISBN: 978-953-6931-79-8.<b><br>