Maj 2015, J. Slavič in L. Knez


```python
# Uvozimo potrebne module

from sympy import *
init_printing()
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```

Vprašanje 1: Na sliki (vir: J. Slavič: Dinamika, meh. nihanja..., 2014) je prikazan trikotnik s stranicami dolžine $a$, $b$, debelino $h$ in gostoto $\rho$.

V simbolni obliki določite masni vztrajnostni moment glede na prikazano os $y$:
$$J_{yy}=\int_0^b y^2\,\rho\,h\,(a-a/b\,y)\,dy.$$
Upoštevajte tudi: $m=a\,b\,h\,\rho/2$. Za izmišljene vrednosti izračunajte numerični rezultat.


```python
a, b, h, rho, y, m = symbols('a, b, h, rho, y, m')
Jyy = (y**2 * rho * h * (a-a/b*y)).integrate((y, 0, b))
Jyy
```

Upoštevanje mase:


```python
Jyy.expand().subs(rho*a*b*h/2, m)
```

Izračunajmo še s številkami:


```python
podatki = {'a': 1, 'b': 2, 'rho': 7800, 'h': 0.01}
```


```python
Jyy.subs(podatki)
```

Vprašanje 2: Izračunajte integral tudi numerično. Uporabite ``scipy.integrate`` in integrirajte glede na pravila: trapezno, Simpsonovo 1/3. Rezultat primerjajte tudi z Gaussovo kvadraturo. Raziščite natančnost in hitrost metod.

Definirajmo najprej funkcjo, ki jo je treba integrirati:


```python
def f(y, p=podatki):
    """
    y: spremenljivka
    p: slovar s podatki: a, b, rho, h
    """
    return y**2 * p['rho']* p['h'] * (p['a']-p['a']/p['b']*y)
```


```python
from scipy import integrate
```

Da dobimo občutek o funkciji, ki jo integriramo, jo prikažimo.


```python
# Definiramo meje integrala 
podatki['ma'] = 0
podatki['mb'] = 2

N = 100
y_d = np.linspace(podatki['ma'], podatki['mb'], N)
f_d = f(y_d, podatki)
plt.plot(y_d, f_d);
```


```python
#%%timeit
dy = y_d[1]-y_d[0]
integrate.trapz(f_d, dx=dy)
```


```python
#%%timeit
integrate.simps(f_d, dx=dy)
```


```python
#%%timeit
integrate.quad(f, 0, podatki['b'], args=podatki)
```

Vprašnje 3: Preštudirajte ``scipy.special.legendre``, ki vam vrne objekt ``orthopoly1d``. Ta objekt ima metodo ``weights``, ki vrne seznam ``[x, w, mu0]`` vrednosti, ki jih uporabimo pri Gaussovi kvadraturi. (Če tukaj vsega ne razumete, ne skrbite preveč, bo asistent pokazal/komentiral). Opazite lahko, da smo vrednosti izpeljali na predavanjih!


```python
from scipy import special
```

Pripravimo Legendrev polinom izbrane stopnje (poskusite več različnih):


```python
poli = special.legendre(3)
```

Pogledamo vozlišča, uteži:


```python
poli.weights
```




    array([[ -7.74596669e-01,   5.55555556e-01,   5.55555556e-01],
           [ -4.78946310e-17,   8.88888889e-01,   8.88888889e-01],
           [  7.74596669e-01,   5.55555556e-01,   5.55555556e-01]])



Vprašanje 4: S pomočjo zgoraj pridobljenih uteži in vozlišč izračunajte integral s pomočjo Gaussove kvadrature: $\sum_iw_i\,f(x_i)$. Pazite na transformacijo mej:
$$x=\frac{b+a}{2}+\frac{b-a}{2}\xi$$
$$\int_a^bf(x)\,dx\approx\frac{b-a}{2}\sum_i A_i\,f(x_i).$$

Pripravimo delni rezultat:


```python
aa = [w*f(((podatki['mb']-podatki['ma'])*x + podatki['mb']+podatki['ma'])/2, podatki) for x, w, mu0  in poli.weights]
aa
```

Upoštevamo transformacijo mej in izračunamo integral:


```python
(podatki['mb']-podatki['ma'])/2*np.sum(aa)
```


```python
integrate.newton_cotes(rn=3)
```




    (array([ 0.375,  1.125,  1.125,  0.375]), -0.0375)




```python

```

81218ff2eb6b9005ac0549e32938690c15c516d2
ipynb
Jupyter Notebook
pypinm-master/vprasanja za razmislek/Vaja 11 - polovica.ipynb
CrtomirJuren/python-delavnica
db96470d2cb1870390545cfbe511552a9ef08720
MIT
pypinm-master/vprasanja za razmislek/Vaja 11 - polovica.ipynb
CrtomirJuren/python-delavnica
db96470d2cb1870390545cfbe511552a9ef08720
MIT
pypinm-master/vprasanja za razmislek/Vaja 11 - polovica.ipynb
CrtomirJuren/python-delavnica
db96470d2cb1870390545cfbe511552a9ef08720
MIT
Qwen/Qwen-72B
1. YES
2. YES