Pochodna (różniczkowanie numeryczne)

Różnica dzielona w przód

Rozwińmy funkcję $f(x)$ w otoczeniu $h$ punktu $x$ w szereg Taylora:

$$ f(x+h) = f(x)+ f'(x)h +O(h^2)$$

Dzieląc przez $h$ otrzymujemy:

$$ \displaystyle \frac{f(x+h)-f(x) }{h}=f'(x) + O(h)$$

Np. różnica centralna:

$$ \displaystyle \frac{f(x)-f(x-h) }{2h}=f'(x) + O(h^2)$$

In [1]:
import numpy as np
x = np.linspace(0, 10, 10)
f = np.sin(x)
f1 = np.cos(x)

df = f[1:] - f[:-1]
dx = x[1:] - x[:-1]

x2 = (x[1:] + x[:-1])/2

fp = df/dx
fp.shape, x.shape


Out[1]:
((9,), (10,))

In [2]:
%matplotlib inline

In [3]:
import matplotlib.pyplot as plt

In [4]:
plt.plot(x2, np.cos(x2), 'o-')
plt.plot(x2, fp, 'ro-')


Out[4]:
[<matplotlib.lines.Line2D at 0x7f428955ac18>]

In [ ]:

Numeryczne całkowanie

Całka to pole, więc obliczmy je w przybliżony sposób za pomocą prostokątów.

Przedział $[a,b]$ jest podzielony na $n$ podprzedziałów $$(a,x_1), (x_1,x_2), ..., (x_{n-1},b),$$ punkt $\xi_{j}$ zaś, dla $j = 1, 2, ..., n$ leży gdzieś wewnątrz $j$-tego odcinka. Tworzymy sumę:

$$S_{n} = \sum_{j=1}^{n}f(\xi_{j})(x_{j} - x_{j-1}) = \sum_{j=1}^{n}f(\xi_{j}) \Delta x_{j}$$

przy czym $x_0 = a$, $x_n = b$ oraz $\Delta x_j = x_j - x_{j-1}$.

Geometrycznie suma ta, nazywana sumą Riemanna, jest sumą pól prostokątów zaznaczonych na rysunku poniżej. Jeżeli będziemy dzielić $[a,b]$ na coraz to więcej odcinków o coraz to mniejszych długościach, otrzymamy jako granicę sum Riemanna całkę Riemanna funkcji $f(x)$, oznaczaną przez:

$$ \int_{a}^{b} f(x)dx = \lim_{|L| \to 0} \sum_{j=1}^{n} f(\xi_{j}) \Delta x_{j} \quad (1) $$

gdzie przez $|L|$ oznaczyliśmy długość największego z podprzedziałów.

Kwadratury

Metody numetyczne obliczania całek sprowadzają sie do obliczania kwadratur - tj. sumy ważonej wartości funkcji w punktach:

$$ \int_a^b f(x)dx = \sum_{i=1}^N w_i f_i,$$

gdzie $f_i = f(x_i)$.

W metodzie trapezów link mamy wagi $w_i=1$ dla wszystkich punków z wyjątkiem: $w_1=w_{N}=\frac{1}{2}$.

Przykład:

$$\int_0^\pi \sin(x)dx = 2$$

In [2]:
import numpy as np 
x = np.linspace(0,np.pi,10)
f = np.sin(x)
w = np.ones_like(x)
w[0] = 0.5
w[-1] = 0.5
h = x[1]-x[0]
print(h*np.sum(w*f))


1.97965081122

Funkcja pierwotna

$$ \int_{a}^{b} f(x)dx = F(b) - F(a) $$$$F(x) = \int_{a}^{x} f(x')dx' $$


In [ ]:

Pochodne funkcji wielu zmiennych


In [1]:
import numpy as np
x = np.linspace(-2,1,40)
y = np.linspace(-2,3,34)
X,Y = np.meshgrid(x,y)

In [2]:
%matplotlib notebook
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

In [ ]:


In [3]:
import sympy
from sympy.abc import x,y
sympy.init_printing(use_latex='mathjax')

In [4]:
from IPython.display import display

In [5]:
f_symb = -sympy.exp(-(x**2+2*y**2))
display(f_symb)

F = sympy.lambdify((x,y),f_symb,np)
Fx = sympy.lambdify((x,y),f_symb.diff(x),np)
Fy = sympy.lambdify((x,y),f_symb.diff(y),np)

fig = plt.figure()
ax = fig.gca(projection='3d')

#ax.plot_wireframe(X,Y,F)
ax.plot_surface(X, Y, F(X,Y), cmap=cm.coolwarm,rstride=1,cstride=1)
#ax.quiver3D(X, Y, Fx, Fy) #, cmap=cm.coolwarm,rstride=1,cstride=1)
ax.set_xlabel('X')
ax.set_ylabel('Y')


$$- e^{- x^{2} - 2 y^{2}}$$
Out[5]:
<matplotlib.text.Text at 0x7fe6ffcc47b8>

In [6]:
x0,y0 = -1.5,-.5
h = 0.2
for i in range(30):
    x0 += -h * Fx(x0,y0)
    y0 += -h * Fy(x0,y0)
    #ax.scatter3D(x0,y0,F(x0,y0),s=240,c='g',marker='o')
    ax.plot([x0,x0],[y0,y0],[-1,F(x0,y0)],c='r')

In [7]:
plt.figure()
plt.contourf(X,Y,F(X,Y))
plt.quiver(X,Y,-Fx(X,Y),-Fy(X,Y))
plt


Out[7]:
<module 'matplotlib.pyplot' from '/opt/conda/lib/python3.5/site-packages/matplotlib/pyplot.py'>

In [8]:
x0,y0 = -1.5,-.5
h = 0.1
for i in range(100):
    x0 += -h * Fx(x0,y0)
    y0 += -h * Fy(x0,y0)
    plt.plot(x0,y0,'go')

hill shading


In [9]:
x_ = np.linspace(-5,5,150)
y_ = np.linspace(-5,5,154)
X,Y = np.meshgrid(x_,y_)

f_symb = sympy.sin(x**2/4+y**2)
expr = sympy.diff(f_symb,x)
display(f_symb)

F = sympy.lambdify((x,y),f_symb,np)
Fx = sympy.lambdify((x,y),f_symb.diff(x),np)
Fy = sympy.lambdify((x,y),f_symb.diff(y),np)


$$\sin{\left (\frac{x^{2}}{4} + y^{2} \right )}$$

In [10]:
plt.figure()
plt.imshow(F(X,Y),alpha=.6,cmap='jet')
plt.imshow(np.diff(F(X,Y),axis=0),alpha=.4,cmap='gray')


Out[10]:
<matplotlib.image.AxesImage at 0x7f2ecc0fb2e8>

In [ ]:


In [ ]:


In [19]:
import sympy
from sympy.abc import x,y
sympy.init_printing(use_latex='mathjax')

In [20]:
from IPython.display import display

In [21]:
f_symb = sympy.sin(x**2+y**2)
lap_f = f_symb.diff(x,2) + f_symb.diff(y,2)
display(lap_f.simplify())


$$- 4 x^{2} \sin{\left (x^{2} + y^{2} \right )} - 4 y^{2} \sin{\left (x^{2} + y^{2} \right )} + 4 \cos{\left (x^{2} + y^{2} \right )}$$

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: