import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import warnings
# Ignore all warnings
warnings.simplefilter('ignore')
Faltung¶
t = np.linspace(-4,5,361)
# funktion f(t)
ft = np.heaviside(t+1,1) - np.heaviside(t-1,1)
# funktion g(t)
def gt(t,a):
return np.heaviside(t-a,1) - np.heaviside(t-a-2,1)
# Darstellung für a = 0
fig, ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(t, ft, label='$f(t)$')
ax.plot(t, gt(t,0), label='$g(t)$', linestyle='--')
ax.set_xlabel('$t$')
ax.legend()
plt.show()
Bestimmung der Faltung $$f(t) \ast g(t) = \int_{-\infty}^\infty f(a) g(t-a) \mathrm{d}a$$
# Achtung: Wechsel der Koordinaten
a = np.linspace(-4,5,361)
fa = np.heaviside(a+1,1) - np.heaviside(a-1,1) # keine Änderung
obs_t = list()
obs_c = list()
for t in [-2,-1/2, 0, +1/2, 1, 1.5, 2, 2.5, 4]:
fig, ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(a, fa, label='$f(a)$')
ax.plot(a, gt(t,a), label='$g(t-a)$', linestyle='--')
leftcorner = max(-1,t-2)
rightcorner = min(1, t)
if leftcorner < rightcorner:
ax.fill_betweenx([0,1], leftcorner, rightcorner, color='green', alpha=0.5)
ax.set_xlabel('$a$')
ax.legend()
ax.set_title(f'Zeitpunkt t = {t}')
ax.grid()
plt.show()
print(f'Fläche der Überschneidung: {max(rightcorner - leftcorner,0)}')
# save observed points
obs_t.append(t)
obs_c.append(max(rightcorner - leftcorner,0))
Fläche der Überschneidung: 0
Fläche der Überschneidung: 0.5
Fläche der Überschneidung: 1
Fläche der Überschneidung: 1.5
Fläche der Überschneidung: 2
Fläche der Überschneidung: 1.5
Fläche der Überschneidung: 1
Fläche der Überschneidung: 0.5
Fläche der Überschneidung: 0
# use numpy convolution
start = -4
end = -start
points_per_unit = 100
t = np.linspace(start,end,(end-start)*points_per_unit+1)
ft = np.heaviside(t+1,1) - np.heaviside(t-1,1)
res = np.convolve(ft, gt(t,0), mode='same')/points_per_unit
fig, ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(t,res)
ax.scatter(obs_t, obs_c, color='red', marker='+', s=150)
ax.grid()
ax.set_xlabel('$t$')
ax.set_ylabel('$f(t) * g(t)$')
plt.show()
$ \begin{align} f(t) \ast g(t) & = \int_{-\infty}^\infty f(a) g(t-a) \mathrm{d}a\\ & = \int_{-\infty}^\infty [\sigma(a+1) - \sigma(a-1)] \cdot [\sigma(t-a) - \sigma(t-a -2)] \mathrm{d}a\\ & = \int_{-\infty}^\infty [\sigma(a+1)\sigma(t-a) - \sigma(a+1)\sigma(t-a -2) - \sigma(a-1)\sigma(t-a) + \sigma(a-1)\sigma(t-a -2) \mathrm{d}a\\ & = \int_{-1}^t \sigma(a+1) \sigma(t-a) \mathrm{d}a - \int_{-1}^{t-2} \sigma(a+1)\sigma(t- a -2) \mathrm{d}a - \int_{1}^t \sigma(a-1)\sigma(t-a) \mathrm{d}a + \int_{1}^{t-2}\sigma(a-1) \sigma(t- a -2) \mathrm{d}a \\ & = (t+1) \sigma(t+1)- (t-1) \sigma(t-1) - (t-1)\sigma(t-1) + (t-3)\sigma(t-3) \\ & = (t+1) \sigma(t+1) - 2(t-1) \sigma(t-1) + (t-3)\sigma(t-3) \end{align} $
start = -4
end = -start
points_per_unit = 100
t = np.linspace(start,end,(end-start)*points_per_unit+1)
ft = np.heaviside(t+1,1) - np.heaviside(t-1,1)
res = np.convolve(ft, gt(t,0), mode='same')/points_per_unit
res_theory = (t+1)*np.heaviside(t+1,1) - 2*(t-1)*np.heaviside(t-1,1) + (t-3)*np.heaviside(t-3,1)
fig, ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(t,res_theory, label='theoretisch')
ax.plot(t,res,label='numerisch', linestyle='--')
ax.grid()
ax.set_xlabel('$t$')
ax.set_ylabel('$f(t) * g(t)$')
ax.legend()
plt.show()
Lösen von DGL¶
Beispiel $$y^\prime + 4y = t^3 \qquad y(1) = 2$$
Lösung $$y(t) = \dfrac{1}{128}\left(32t^3 - 24t^2 + 12t - 3\right) + 101.94 \mathrm{e}^{-4t}$$
t = np.linspace(0,1.5,201)
y = 1/128*(32*t**3 - 24*t**2 + 12*t - 3) + 101.94*np.exp(-4*t)
y_prime = t**3 - 4*y
y_prime2 = 1/128*(3*32*t**2 - 2*24*t**1 + 12) + (-4)*101.94*np.exp(-4*t)
fig, axs = plt.subplots(1,2,figsize=(12,4))
cax = axs[0]
cax.plot(t,y)
cax.scatter([1], [2], marker='x', s=100, color='red')
cax.set_xlabel('$t$')
cax.set_ylabel('$y(t)$')
cax.grid()
cax = axs[1]
cax.plot(t, y_prime, label='aus DGL')
cax.plot(t, y_prime2, linestyle='--', label='Ableitung der Lsg')
cax.set_xlabel('$t$')
cax.set_ylabel('$y^\prime(t)$')
cax.legend()
cax.grid()
plt.show()