4.6. Temperaturgleichung#

Die Lösung einer gewöhnlichen Differentialgleichung, hier der Temperaturentwicklung, wird in dieser Aufgabe nochmals vorgestellt. Dabei wird die Gleichung erweitert und ein komplexer Quellterm implementiert.

Aufgabenteil A#

In der Vorlesung wurde folgende Differentialgleichung vorgestellt.

\[ \sf T'(t) = -(T-T_\infty) + Q(t) \]

Verwenden Sie die Funktion scipy.integrate.solve_ivp um die Gleichung für \(\sf T_0 = 20\) mit \(\sf T_\infty=10\) und \(\sf Q=0\) zu lösen. Die Lösung soll bis \(\sf t=25\) berechnet werden. Erstellen und beschriften Sie den Lösungsplot.

Lösungshinweis#

Die Ausgabe könnte wie folgt aussehen.

Lösungsvorschlag#

Hide code cell content
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
Hide code cell content
# Rechte-Hand-Seite der Gleichung, hier mit Q=0
def f_T1(t, T):
    Tinfty = 10
    return -(T-Tinfty)
Hide code cell content
# Anfangswert der Lösung
T0 = 20

# Aufruf des Lösers
res1 = scipy.integrate.solve_ivp(f_T1, [0, 25], [T0])
Hide code cell content
plt.plot(res1.t, res1.y[0], marker='.')
plt.xlabel("Zeit t")
plt.ylabel("Temperatur T")
plt.grid()

# Ausgabe für den Lösungshinweis
# plt.savefig('teil1.png')
../../../../_images/58d9a99469d91fc95fe2ba2da82c71856d3521b36e2e631a168e86cbc2090ba2.png

Aufgabenteil B#

Fügen Sie einen neuen Koeffizienten \(\sf \beta\) als Vorfaktor vor den ersten Term auf der rechten Seite der Gleichung. Lösen Sie nun die Gleichung für \(\sf \beta=0.2\). Stellen Sie graphisch die Lösungen mit und ohne das \(\sf \beta\) in einer Graphik dar.

Lösungshinweis#

Die Ausgabe könnte wie folgt aussehen.

Lösungsvorschlag#

Hide code cell content
def f_T2(t, T):
    Tinfty = 10
    beta = 0.2
    return -beta*(T-Tinfty)
Hide code cell content
T0 = 20
res2 = scipy.integrate.solve_ivp(f_T2, [0, 25], [T0])
Hide code cell content
plt.plot(res1.t, res1.y[0], marker='.', label='beta=1')
plt.plot(res2.t, res2.y[0], marker='.', label='beta=0.2')

plt.xlabel("Zeit t")
plt.ylabel("Temperatur T")
plt.grid()
plt.legend()

# Ausgabe für den Lösungshinweis
# plt.savefig('teil2.png')
<matplotlib.legend.Legend at 0x7f9171639660>
../../../../_images/e12b5103f17d43d32b5d3d5d57f0e997a984b3b2fd6d524855e28b329202aebf.png

Aufgabenteil C#

Nun wird der Quellterm \(\sf Q(t)\) erweitert. Dieser soll für die Zeit zwischen 5 und 25 zwei volle Sinus-Schwingungen mit einer Amplitude von 10 durchführen, sonst ist er Null. Lösen Sie die Gleichung mit \(\sf \beta=0.2\) bis zur Zeit \(\sf t=50\). Stellen Sie neben der Lösung auch die Quelltermfunktion \(\sf Q(t)\) in einer gemeinsamen Graphik dar.

Lösungshinweis#

Die Ausgabe könnte wie folgt aussehen.

Lösungsvorschlag#

Hide code cell content
def f_Q(t):
    Q = 0
    if t >= 5 and t <= 25: 
        Q = 10 * np.sin((t-5)/20 * 4 * np.pi)
    return Q
Hide code cell content
def f_T3(t, T):
    Tinfty = 10
    return -0.2*(T-Tinfty) + f_Q(t)
Hide code cell content
T0 = 20
res3 = scipy.integrate.solve_ivp(f_T3, [0, 50], [T0])
Hide code cell content
plt.plot(res3.t, res3.y[0], marker='.', label='T(t)')

Q = np.zeros_like(res3.t)
for i in range(len(Q)):
    Q[i] = f_Q(res3.t[i])
    
plt.plot(res3.t, Q, marker='.', label='Q(t)')

plt.xlabel("Zeit t")
plt.ylabel("Funktionswert")
plt.grid()
plt.legend()

# Ausgabe für den Lösungshinweis
# plt.savefig('teil3.png')
<matplotlib.legend.Legend at 0x7f916f48b040>
../../../../_images/53104b6931c575c4baf2f6abe0dfb45e54902ab86f2f9d98241011445ace88bf.png

Aufgabenteil D#

Die graphische Darstellung, siehe Lösungshinweis, deutet auf eine sehr sprunghafte Lösung. Verwenden Sie das optionale Argument max_step der Funktion scipy.integrate.solve_ivp um einen kleineren Zeitschritt, z.B. \(\sf \Delta t=0.25\), zu wählen. Vergleichen Sie nun die beiden Lösungen.

Lösungshinweis#

Die Ausgabe könnte wie folgt aussehen.

Lösungsvorschlag#

Hide code cell content
T0 = 20
res4 = scipy.integrate.solve_ivp(f_T3, [0, 50], [T0], max_step=0.5)
Hide code cell content
plt.plot(res3.t, res3.y[0], marker='.', label='T(t) – dt variabel')
plt.plot(res4.t, res4.y[0], marker='.', label='T(t) – dt max 0.5')

Q = np.zeros_like(res4.t)
for i in range(len(Q)):
    Q[i] = f_Q(res4.t[i])
    
plt.plot(res4.t, Q, marker='.', label='Q(t)')

plt.xlabel("Zeit t")
plt.ylabel("Funktionswert")
plt.grid()
plt.legend()

# Ausgabe für den Lösungshinweis
# plt.savefig('teil4.png')
<matplotlib.legend.Legend at 0x7f91718b7cd0>
../../../../_images/05758a7cc008640ec5aa74ad6be7db973a2756bdb9343b33deb01e07f76b7ffc.png