Numerik von Differenzialgleichungen - home.hs …weth0002/veranstaltungen/details/docs/... ·...

Thomas

Westermann

ComputersimulationMit einer Einfuhrung in die Systemanalyse:

Modellbildung, Simulation, Interpretation

c© Thomas Westermann 20.4.2017

Thomas WestermannHochschule Karlsruhe - Technik und WirtschaftMoltkestr. 30Fakultat EIT76133 Karlsruhe

Homepage zur Veranstaltung:

http://www.home.hs-karlsruhe.de/˜weth0002/ —> Veranstaltungen —> Computersimulation

Email: [email protected]

c© Das Copyright samtlicher Materialien liegt beim Autor. Materialien durfen zum privaten Gebrauchoder zu Studienzwecken heruntergeleden werden. Eine Vervielfaltigung, Verbreitung oder Bereitstellungfur Dritte bedarf der ausdrucklichen Genehmigung des Autors.

Inhaltsverzeichnis

1 Simulation des Frequenzverhaltens fur RCL-Schaltungen . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Berechnung der Ubertragungsfunktion UA

U0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.1.1 Spannungsteiler = Kette aus einem Glied . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.1.2 Kette aus zwei Gliedern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.1.3 Kette aus drei Gliedern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.1.4 Beispiele mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Dimensionierung von Hoch- und Tiefpassen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.2.1 Beispiele mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.3 Ortskurven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.4 Bode-Diagramm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2 Simulation des Zeitverhaltens fur elektrische RCL-Schaltungen . . . . . . . . . . . . . . . . . . . . 172.1 Physikalische Gesetzmaßigkeiten der Bauelemente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Aufstellen der DG einer Schaltung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 Losen der DG mit dem Euler-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.4 Einfache Beispiele mit L, C,R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.1 Reihenschaltung von Spule L und Widerstand R . . . . . . . . . . . . . . . . . . . . . . . . . . 192.4.2 Reihenschaltung von Kapazitat C und Widerstand R . . . . . . . . . . . . . . . . . . . . . . . 192.4.3 Reihenschaltung von R, C und L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.4.4 Parallelschwingkreis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.4.5 Parallelschwingkreis mit realen Spulen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.5 Aufstellen der DG fur Filterschaltungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.5.1 Aufstellen der DG fur einen Tiefpass (TP2PiCLC) . . . . . . . . . . . . . . . . . . . . . . . . . 222.5.2 Beispiele mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.6 Losen der DG fur den Tiefpass TP2PiCLC mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3 Simulation von mechanischen Systemen mit Massenpunkten . . . . . . . . . . . . . . . . . . . . . . 29

4 Simulation geladener Teilchen in elektromagnetischen Feldern . . . . . . . . . . . . . . . . . . . . . 33

5 Signal- und Systemanalyse mit der FFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.1 Grundlegende Formeln . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.2 Eigenschaften der DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.3 Berechnung der DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 395.4 FFT mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.4.1 Signalanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.4.2 Systemanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

VI Inhaltsverzeichnis

6 Signal- und Systemanalyse mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.1 FFT eines Sinussignals mit Maple: Wie hangt der Spektralbereich von der Anzahl der

Abtastpunkte ab? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.2 Das Hanning-Fenster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 456.3 Die Prozedur plot fft mit der Option Hanning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476.4 Maximal erfassbare Frequenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496.5 Frequenzauflosung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 506.6 Analyse von Messdaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 516.7 System-Analyse mit Maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

1. Simulation des Frequenzverhaltens furRCL-Schaltungen

Im Folgenden werden wir eine Systematik entwickeln, wie man die komplexe Ubertragungsfunktionfur beliebig lange Ketten sukzessive berechnet, wenn die Schaltung aus R−, C− oder L−Elementenaufgebaut ist. Im Anschluss daran werden einige Spezialfalle fur Tiefpass, Hochpass, Bandpass undBandsperren diskutiert. Das Ziel ist es anschließend, nicht nur zu gegebener Schaltung die Ubertra-gungsfunktion zu bestimmen, sondern zu gegebener Grenzfrequenz fur einen Hoch- oder Tiefpass eineDimensionierung der einzelnen Elemente der Kette vorzunehmen. Dabei soll die UbertragungsfunktionUA

U0einen moglichst glatten Verlauf im Arbeitsbereich besitzen.

Qualitativ werden fur die Filter folgende Verlaufe von | UA

U0| gefordert:

Abb. 1.1. Ubertragungsverhalten von Tiefpass, Hochpass, Bandpass und Bandsperre

Zur Realisierung werden wir nur Ketten betrachten, die sich aus Π- oder T -Gliedern zusammensetzen:

Abb. 1.2. Schaltungen, die auf T-Gliedern basieren

Abb. 1.3. Schaltungen, die auf Π-Gliedern basieren

2 1 Simulation des Frequenzverhaltens fur RCL-Schaltungen

Die Ubertragungsfunktion ist definiert als das Verhaltnis von Ausgangs- zu Eingangsspannung, wenndie Eingangsspannung als komplexe Wechselspannung U0 = U0 eiωt angenommen wird:

H(ω) =UA

U0

.

Diese komplexe Ubertragungsfunktion lasst sich dann graphisch uber den Betrag |H(ω)| und die Phaseϕ(ω) darstellen:

H(ω) = |H(ω)|eiϕ(ω)

mit |H(ω)| =√

H(ω) ·H∗(ω) und tan(ϕ(ω)) = Im(H(ω))Re(H(ω)) .

Man beachte, dass fur den folgenden Algorithmus die einzelnen Elemente bzw. Kettenglieder noch nichtspezifiziert werden; d.h. die Formeln gelten sowohl fur Π- als auch fur T -Glieder gleichermaßen. Furdie Impedanzen (Quer- und Langsimpedanzen) muss dann konkret fur die jeweilige Schaltung noch

RΩ (ohmscher Widerstand)iωL (Impedanz einer Spule mit Induktivitat L)

1iωC (Impedanz eines Kondensators mit Kapazitat C)

gesetzt werden.

Tabelle: Filterschaltungen

Tiefpasse aus T -Glieder bestehend:

Hochpasse aus T -Glieder bestehend:

Tiefpasse aus Π-Glieder bestehend:

1 Simulation des Frequenzverhaltens fur RCL-Schaltungen 3

Hochpasse aus Π-Glieder bestehend:

Bandpasse aus Π-Glieder bestehend:

Bandpasse aus Π-Glieder bestehend:

Bandpasse aus T -Glieder bestehend:

Bandsperren aus T -Glieder bestehend:



U0

1.1.1 Spannungsteiler = Kette aus einem Glied

z1 = Langsimpedanzy1 = Querimpedanz

UA

U0

=y1I

z1I + y1I=

y1

z1 + y1

1.1.2 Kette aus zwei Gliedern

Wir betrachten zunachst die zweite Masche und ersetzendiese durch einen Ersatzwiderstand

y1p =1

1y1

+ 1z2+y2

=y1(z2 + y2)y1 + z2 + y2

Fur die verbleibende Masche konnen wir die Spannungstei-lerregel anwenden:

U1

U0

=y1p

y1p + z1

Damit kennen wir die Spannung, die bei y1p abfallt undkonnen nun nochmals die Spannungsteilerregel auf die 2.Masche anwenden:

U2

U1

=y2

y2 + z2

Zusammenfassung:z1, z2 Langsimpedanzeny1, y2 Querimpedanzen

y1p = y1 · (y2 + z2)/(y1 + y2 + z2)u1 = u0 · y1p/(z1 + y1p)u2 = u1 · y2/(z2 + y2)


U05

1.1.3 Kette aus drei Gliedern

Man ersetzt die Schaltung sukzessive von rechtsnach links durch Ersatzwiderstande, so dass zumSchluss nur noch eine Masche mit einem Paral-lelwiderstand ubrig bleibt.

y2p =y2(y3 + z3)y2 + y3 + z3

y1p =y1(y2p + z2)y2 + y2p + z2

Anschließend wendet man nun von links nachrechts die Spannungsteilerregel an

→ →

U1

U0

=y1p

y1p + z1→ U2

U1

=y2p

y2p + z2→ U3

U2

=y3

y3 + z3.

Zusammenfassungz1, z2, z3 Langsimpedanzeny1, y2, y3 Querimpedanzen

y2p = y2 · (y3 + z3)/(y2 + y3 + z3)y1p = y1 · (y2p + z2)/(y1 + y2p + z2)U1 = U0 · y1p/(y1p + z1)U2 = U1 · y2p/(y2p + z2)U3 = U2 · y3/(y3 + z3)


1.1.4 Beispiele mit Maple

Wir berechnen die komplexe Ubertragungsfunktion fur (1) einen Hochpass, (2) einen Tiefpass, (3)einen Bandpass und (4) eine Bandsperre mit Maple. Die zugehorigen Maple-Worksheets konnen durchAnklicken gestartet und anschließend ausgefuhrt werden.

HP2TCLC Hochpass, der aus zwei T -Gliedern besteht.TP2PiCLC Tiefpass, der aus zwei Π-Gliedern besteht.BP2PiLCp Bandpass, der aus zwei Π-Gliedern besteht.BS2TLCp Bandsperre, die aus zwei T -Gliedern besteht.


U07

(1.) Hochpass: 2-T-CLC-Glieder HP2TCLC

-------R-----C-------------1/2 C-----------C-------| | |L L R| | |

---------------------------------------------------

> Z1:= R + 1/(I*w*C):

> Z2:= 1/(I*w*C/2):

> Z3:= 1/(I*w*C):

> Y1:= I*w*L:

> Y2:=I*w*L:

> Y3:=R:

> Y2p:= Y2*(Y3+Z3)/(Y2+Y3+Z3):

> Y1p:= Y1*(Y2p+Z2)/(Y1+Y2p+Z2):

> u1:= Y1p/(Y1p+Z1):

> u2:= u1*Y2p/(Y2p+Z2):

> u3:= u2*Y3/(Y3+Z3):

> wiu3:=argument(u3):

> beu3:=abs(simplify(u3));

beu3 :=12

∣∣w5 L2 C3 R/(I w5 L2 C3 R + 2 w4 L2 C2 − 4 I w3 LC2 R− 3 w2 LC + R2 w4 C3 L

−R2 w2 C2 + 2 I R w C + 1)∣∣

Die Ubertragungsfunktion ist eine gebrochenrationale Funktion in w mit hochstem Exponent 5, da dieSchaltung 5 Energiespeicher besitzt. Die graphische Darstellung des Betrags der Ubertragungsfunktionerhalt man durch> R:=1000: C:=5.28e-9: L:=3.128e-3:

> titel:=‘Ubertragungsfunktion HP2TLCL‘:

> plot( beu3,w=0..400000,title=titel,thickness=2);

bertragungsfunktion HP2TLCL

0

0.1

0.2

0.3

0.4

0.5

100000 200000 300000 400000w

Fur R=1000 Ohm, C=5.28e-9 F, L=3.128e-3 H erhalt man die Grenzfrequenz wg=175000 1/s.Die graphische Darstellung der Phase erhalt man durch> plot( wiu3,w=0..400000,thickness=2);


(2.) Tiefpass: 2-Pi-CLC-Glieder: TP2PiCLC

-------R----------------------- L-----------------L-------------| | | |C 2 C C R| | | |

----------------------------------------------------------------

(kr=0.8; kf=1.3)

> Z1:= R:

> Z2:= I*w*L:

> Z3:= I*w*L:

> Y1:= 1/(I*w*C):

> Y2:= 1/(I*w*2*C):

> Y3:= 1/(I*w*C+1/R):

> Y2p:= Y2*(Y3+Z3)/(Y2+Y3+Z3):

> Y1p:= Y1*(Y2p+Z2)/(Y1+Y2p+Z2):

> u1:= Y1p/(Y1p+Z1):

> u2:= u1*Y2p/(Y2p+Z2):

> u3:= u2*Y3/(Y3+Z3):

> wiu3:=argument(u3):

> beu3:=abs(simplify(u3));

beu3 :=12

∣∣R/(−I R + 4 I w2 LC R + w L− 2 I R w4 L2 C2 − w3 L2 C + 2 w C R2

− 3 R2 C2 Lw3 + w5 L2 C3 R2)∣∣

> R:=2500: C:=1e-9: L:=10e-3:

> titel:=‘Ubertragungsfunktion TP2PiCLC‘:

> plot( beu3,w=0..1000000,title=titel,thickness=2);

> titel:=‘Phase TP2PiCLC‘:

> plot( wiu3,w=0..1000000,title=titel,thickness=2);

bertragungsfunktion TP2PiCLC

0

0.1

0.2

0.3

0.4

0.5

200000 400000 600000 800000 1e+06w

Phase TP2PiCLC

–3

–2

–1

0

1

2

3

200000 400000 600000 800000 1e+06w

Die Grenzfrequenz ist wg=445000 1/s.

1.2 Dimensionierung von Hoch- und Tiefpassen 9

1.2 Dimensionierung von Hoch- und Tiefpassen

Bisher haben wir diskutiert, wie man zu einer gegebenen Kette die Ubertragungsfunktion berechnetund dabei die Grenzfrequenz ωg bestimmt. Jetzt werden wir den anspruchsvolleren Weg: Wir geben dieGrenzfrequenz ωg und den Widerstand R vor. Gesucht sind dann die passenden Werte fur L und C.

Wenn man die folgenden Voraussetzungen macht, kann man durch gezieltes Probieren mit Maple sehrgute Resultate gewinnen.

Voraussetzungen:(1) Der Aufbau der Kette erfolgt entweder durch identische Π− oder T−Glieder.

(2) Alle Y - und Z- Impedanzen bestehen aus reinen Blindwerten (also nur L und C; kein R). Diesist zwar nicht realisierbar, aber es kann zunachst von verlustlosen Spulen ausgegangen werden.Nachdem angepasste Betriebsparameter gefunden sind, kann der Einfluss von Spulenwiderstandenberucksichtigt werden.

(3) Der Innenwiderstand Ri ist gleich dem Lastwiderstand RA (Leistungsanpassung der Schaltung bzw.Spannungsteilerregel).

Generelle Formeln

(1) Skalierung der Grenzfrequenz: Die Grenzfrequenz ωg bei Hoch- oder Tiefpass ist definiert als dieFrequenz bei der |H(ω)| auf 1√

2des Durchlasswertes abnimmt. Es gilt der Zusammenhang

ωg = kf√LC

, (1.1)

wenn kf eine schaltungsspezifische Konstante darstellt.

(2) Skalierung des Widerstandes: Der Zusammenhang zwischen Widerstand und L, C lautet

R = kR ·√

LC , (1.2)

wenn kR eine schaltungsspezifische Konstante darstellt.

Um die schaltungsspezifischen Konstanten kf und kR zu bestimmen, wahlen wir fur eine ersteSimulation L = C = 1. Denn dann entnimmt man bei gegebener Schaltung den angepasstenWiderstand R bei dem die Ubertragungsfunktion einen moglichst glatten Verlauf besitzt. Damiterhalt man den Skalierungsfaktor kR. Aus dem Ubertragungsverhalten liest man den Skalierungs-faktor kf ab, der im Falle L = C = 1 genau der Grenzfrequenz entspricht.


Vorgehensweise zur Dimensionierung:

(1) Simulation Man setze L = C = 1 und variiere R1 = RA solange bis ein moglichst glatterFrequenzgang gefunden wird ⇒ Ran (angepasster Widerstand). Gesucht sind dann zu gegebenerKreisfrequenz ωg und vorhandenem Widerstand R die zugehorigen Großen von L und C. Die Uber-tragungsfunktion wird nun so skaliert werden, dass das Ubertragungsverhalten gleich bleibt:

(2) Skalierung von L und C bei gegebenem ωg und R

(1.1) · (1.2) ωg ·R = kR · kf ·1C⇒ C = kf ·kR

ωg·R (1.3)

(1.2)(1.1)

R

ωg=

kR

kf· L ⇒ L = kf ·R

ωg·kR(1.4)

(3) Bestimmt man zu den berechneten Werten wiederum die Ubertragungsfunktion, dann hat die qua-litativ den gleichen Verlauf wie fur die Parameter L = C = 1, nur jetzt mit der Grenzfrequenzωg.

Beispiel: Gesucht ist ein Hochpass bestehend aus einem Π-Glied, der bei einem Widerstand von R =500Ω eine Grenzfrequenz von ωg = 1000s−1 besitzt:

Abb. 1.4. Π-Glied

Mit

z1 = R, z2 =1

iωC, y1 = iωL, y2 =

11R + 1

iωL

folgt

y1p =y1(y2 + z2)y1 + y2 + z2

undu1 =

y1p

y1p + z1

u2 = u1 ·y2

y2 + z2.

(1) Um die schaltungsspezifischen Konstanten kf und kR zu bestimmen, setzen wir L = C = 1. DurchVariation von R im Bereich R = 0.5 . . . 1.5 erhalten wir fur R = 0.8 einen glatten Verlauf derUbertragungsfunktion. Aus dem Schaubild der Ubertragungsfunktion lesen wir die Grenzfrequenzωg = 0.57 ab.

⇒ kR = 0.8 und kf = 0.57


(2) Damit ergeben sich die Dimensionierungen fur L und C:

(I) : C =0.57 · 0.81000 · 500

= 0.912 · 10−6F

(II) : L =0.57 · 5000.8 · 1000

= 0.355H


In den beiden folgenden Worksheets werden einen Hochpass und ein Tiefpass skaliert. Die zugehorigenMaple-Worksheets konnen durch Anklicken gestartet und anschließend ausgefuhrt werden.

(1.) Tiefpass: Gegeben sei ein Tiefpass, der aus zwei Π-Glieder zusammengesetzt ist. Gesucht sind Lund C, so dass er bei R = 1000Ω eine Grenzfrequenz von ω = 200001/s besitzt.

(2.) Hochpass: Gegeben ist ein Hochpass, der aus einem Π-Glied zusammengesetzt ist. Gesucht sind Lund C, so dass er bei R = 500Ω eine Grenzfrequenz von ω = 10001/s besitzt.


Beispiel mit Maple: Skalierung eines Tiefpasses.

Gegeben sei ein Tiefpass, der aus zwei Pi-Glieder zusammengesetzt ist. Um die Ubertragungsfunktion desTiefpasses zu berechnen, verwenden wir die Prozedur kette. Hierzu setzen wir fur die Langsimpedanzen> Z[1] := R: Z[2] := I*w*L: Z[3] := I*w*L:

und fur die Querimpedanzen> Y[1] := 1/(I*w*C): Y[2] := 1/(I*w*2*C): Y[3]:= 1/(I*w*C +1/R ):

Als Eingangsspannung wahlen wir> U0 := 1:

Da der Tiefpass aus 3 Gliedern besteht, ist n=3 und der Aufruf lautet> kette(U0,Z,Y,3);

−12

I R/(−I R + 4 I w2 LR C + w L− 2 I R w4 L2 C2 − w3 L2 C + 2 R2 w C − 3 w3 LC2 R2

+ w5 L2 C3 R2)

Die Ubertragungsfunktion ist eine komplexe gebrochenrationale Funktion in w mit dem hochsten auf-tretenden Exponent 5: Die Kette enthalt 5 unabhangige Energiespeicher. Zur Dimensionierung der L-und C-Teile setzen wir zunachst> C:=1: L:=1:

und zeichnen die Ubertragungsfunktion 3-dimensional> plot3d(abs(H), w=0..2, R=0..2,axes=boxed);

00.5

11.5

2

w

00.5

11.5

2

R

00.10.20.30.40.5


Die Ubertragungsfunktion hat bei R ≈ 0.9 einen glatten Verlauf. Wir zeichnen daher H(w) fur 0.7 ≤R ≤ 1.1 detaillierter:

> with(plots):

> r:=0.7: dR:=0.1:

> for i from 1 to 4

> do

> titel:=convert(r,string):

> p||i:=plot( abs(subs(R=r,H)), w=0..2, thickness=2,title=titel);

> r:=r + dR:

> od:

> display([seq(p||i,i=1..4)],insequence=true);

Aus dem Graphen entnimmt man die Grenzfrequenz> kf := 1.4 :

bei dem optimalen Widerstand> kr := 0.8 :

Bestimmung von L und C fur R=500 und w g=20000:> R := 1000 : wg := 20000 :

> C := kf ∗ kr/(wg ∗ R);

> L := kf/kr ∗ R/wg;

C := .5600000000 10−7 L := .08750000000

Die Ubertragungsfunktion fur diese angepassten Dimensionierungen ist:> plot( abs(H), w=0..40000, thickness=2, title=‘Ubertragungsfunktion‘);

bertragungsfunktion

0

0.1

0.2

0.3

0.4

0.5

10000 20000 30000 40000w


1.3 Ortskurven

Eine zu Betrag und Phase alternative graphische Darstellung der Ubertragungsfunktion ist die sog.Ortskurve. Bei dieser Darstellung von H(ω) werden der Betrag |H(ω)| und der Phasenwinkel ϕ(ω)gleichzeitig in ein Schaubild im Polar-Koordinatensystem gezeichnet. Der Betrag entspricht dem Radiusund der Phasenwinkel dem Winkel. Dabei variieren |H(ω)| und ϕ(ω) gleichzeitig als Funktion von ω.Fur unser Beispiel des Hochpasses mit einem Π-Glied ist

0.1

0.2

0.3

0.4

–0.2 –0.1 0.1 0.2 0.3

(a) L = C = 1. ω variiert von 0 . . . 2.

0

0.1

0.2

0.3

0.4

–0.2 –0.1 0.1 0.2 0.3

(b) L = 0.4H, C = 0.9 · 10−6F . ω variiertvon 0 . . . 3000.

Abb. 1.5. Ortskurven

Die beiden Ortskurven sind fur die entsprechenden ω-Parameterbereiche identisch. In Maple wird dieOrtskurve dargestellt entweder durch den Befehl> plot([Re(H(ω)), Im(H(ω)), ω = 0 . . . 2]);

wobei Re(H(ω)) bzw. Im(H(ω)) fur den Realteil bzw. Imaginarteil der komplexen Funktion H(ω)steht oder> with(plots):

> polarplot([abs(H(ω)), argument(H(ω)), ω = 0..2]);

wenn abs(H(ω)) der Betrag und argument(H(ω)) die Phase von H(ω) bedeutet.

1.3 Ortskurven 15

Fur die anfangs des Kapitels diskutierten Beispiele des (1) Hoch-, (2) Tief-, (3) Bandpasses bzw. (4)der Bandsperre ergeben sich die folgenden Ortskurven:

–0.4

–0.2

0

0.2

0.4

–0.4 –0.2 0.2

(a) Hochpass HP2TLCL

–0.4

–0.2

0.2

0.4

–0.4 –0.2 0.2 0.4

(b) Tiefpass TP2PiCLC

–0.4

–0.2

0

0.2

0.4

–0.4 –0.2 0.2 0.4

(c) Bandpass BP2PiLCr

–0.4

–0.2

0

0.2

0.4

–0.4 –0.2 0.2 0.4

(d) Bandsperre BS2TLCp

Abb. 1.6. Ortskurven

Beobachtung des Phasenwinkels: Beim Variieren der Frequenz von 0 bis∞ wandert der Phasenwinkeluber N · 90. Dabei ist N die Anzahl der Energiespeicher.


1.4 Bode-Diagramm

Eine weitere Moglichkeit das Ubertragungsverhalten H(ω) graphisch darzustellen besteht darin denBetrag der Ubertragungsfunktion, |H(ω)|, oder die Phase, argument(H(ω)), als Funktion von ω zuzeichnen.> plot(abs(H(ω)), ω = 0..ωmax);

> plot(argument(H(ω)), ω = 0..ωmax);

Zu dieser Darstellung hat sich in den Anwendungen eine alternative Moglichkeit durchgesetzt, namlichder 10-er Logarithmen des Betrags bei logarithmischer Skalierung der ω-Achse zu zeichnen, dem sog.Bode-Diagramm. Genau gesagt wird

20 dB log10(|H(ω)|)

bei logarithmischer Skalierung der ω-Achse aufgetragen. In Maple wird das Bode-Diagramm realisiertdurch den Befehl> with(plots);

> semilogplot(20*log[10](abs(H(ω)), ω = 0.1..ωmax);

Man beachte dabei, dass ω nicht bei 0 beginnen darf, da eine logarithmische Skalierung der ω-Achsevorgenommen wird und log(0) nicht definiert ist!

6.35

6.30

6.25

6.20

6.15

6.10

6.05

1. .1e2 .1e3 .1e4 .1e5 1e+05 1e+06w

Abb. 1.7. Bode-Diagramm

2. Simulation des Zeitverhaltens fur elektrischeRCL-Schaltungen

Im Folgenden werden fur elektrische Schaltungen, die sich aus RCL-Gliedern zusammensetzen die zu-gehorigen Differenzialgleichungen aufgestellt, die dann mit Hilfe des Rechners gelost werden. Das Zielist fur beliebige Eingangssignale die zugehorigen Ausgangssignale zu berechnen.

2.1 Physikalische Gesetzmaßigkeiten der Bauelemente

Widerstand: U(t) = R · I(t)

Der Strom durch den Widerstand R ist I(t) = 1RU(t), wenn U(t) der Spannungsabfall am Widerstand

R.

Spule mit Induktivitat L: U(t) = LdI(t)dt

Fließt durch eine Spule L ein Strom I, so ist der Spannungsabfall an der Spule U proportional dIdt . Die

Proportionalitatskonstante bezeichnet man mit Induktivitat L: U(t) = LdI(t)dt .

Kondensator mit Kapazitat C: I(t) = C dU(t)dt

Liegt am Kondensator die Spannung U(t), so ist die auf dem Kondensator liegende Ladung Q(t) =C · U(t), wenn C die Kapazitat des Kondensators ist. Wegen I = dQ

dt folgt fur den Strom durch den

Kondensator I(t) = C · dU(t)dt . Naturlich fließt der Strom nicht durch den Kondensator, sondern auf der

einen Seite fließt Ladung zu, auf der anderen Seite fließt Ladung ab.

2.2 Aufstellen der DG einer Schaltung

(1) Jedem Energiespeicher wird eine Zustandsvariable zugeordnet:Jedem Kondensator Ci wird seine Spannung Ui, jeder Spule Li wird ihr Strom Ii alsZustandsvariable zugeordnet. Den Widerstanden wird keine Zustandsvariable zugeordnet, daWiderstande keine Energiespeicher, sondern nur Energieverbraucher darstellen.

(2) Der Maschensatz und der Knotensatz werden auf die Schaltung angewendet. Das Ziel ist fur je-de Zustandsvariable eine DG 1. Ordnung zu erhalten. Sind Spule und Kondensator in Reihegeschaltet, wird formal ein weiterer Knoten eingefuhrt.

(3) In der Regel fuhrt die Vorgehensweise zu je einer DG fur jede Zustandsvariable, die außer der

Ableitung der Zustandsvariablen keine weiteren Ableitungen enthalt. 4! Achtung: Manchmalsind aber mehrere Ableitungen 1. Ordnung in einer DG vorhanden. Dann muss noch ein linearesGleichungssystem gelost werden, damit man die gewunschte Struktur erhalt.

18 2 Simulation des Zeitverhaltens fur elektrische RCL-Schaltungen

2.3 Losen der DG mit dem Euler-Verfahren

Gegeben ist die DG mit Anfangsbedingung

y′(t) = f(t, y(t)) + Ue(t)y(0) = y0.

Nach dem Eulerschen-Streckenzug-Verfahren kann diese DG schrittweise gelost werden:

y(ti+1) = y(ti) + y′(ti) · dt

bzw. kurzyneu = yalt + (y′)alt · dt.

Abb. 2.1. Euler-Verfahren

Dies liefert bei gegebener Funktion Ue(t) und Anfangsbedingung y(0) = y0 den folgenden Algorithmuszum numerischen Losen der DG:

y = y0for t from tmin by dt to tmax do

y = y + (f(t,y) + Ue(t)) * dtend do

Um anschließend die numerische Losung zeichnen zu konnen, speichern wir zu jedem Zeitpunkt sowohlden Zeitpunkt t als auch den Funktionswert y in einer Liste ab. Der so erweiterte Algorithmus lautetdann

y=y0data[0] = [tmin, y]i=1

for t from tmin by dt to tmax doy = y + (f(t,y) + Ue(t)) * dtdata[i] = [t+dt, y]i = i + 1

end do

2.4 Einfache Beispiele mit L, C, R 19

2.4 Einfache Beispiele mit L, C, R

2.4.1 Reihenschaltung von Spule L und Widerstand R

Energiespeicher: LZustandsvariable: I

Maschensatz: U0 = UR + UL = R · I + LdIdt

⇒ dIdt = (U0 −RI)/L → DG fur I

Algorithmus: U0 = A · sin(2πft)i = 1I = 0for t from tmin by dt to tmax

doI = I + (U0 −RI)/L · dtdata[i] = [t + dt, I]i = i + 1

end do

2.4.2 Reihenschaltung von Kapazitat C und Widerstand R

Energiespeicher: CZustandsvariable: U

Maschensatz: U0 = UR + U = R · I + U = R · C dUdt + U

⇒ dUdt = (U0 − U)/(R · C) → DG fur U

Algorithmus: U0 = A · sin(2πft)U = 0for t from tmin by dt to tmax

doU = U + (U0 − U)/(R · C) · dt...

end do


2.4.3 Reihenschaltung von R, C und L

C

R L

I

U

K

U0

Zwei Energiespeicher: L,CZwei Zustandsvariable: I, U

Maschensatz: U0 = UR + UL + U = R · I + LdIdt + U

Knotensatz fur K: I = C · dUdt

⇒dIdt = (U0 − U −RI)/L → DG fur I

dUdt = I/C → DG fur U

Algorithmus: U0 = A · sin(2πft)U = 0 und I = 0for t from tmin by dt to tmax

doI = I + (U0 − U −RI)/L · dtU = U + I/C · dt...

end do

2.4.4 Parallelschwingkreis

C

R

L

I

U

K

U0

M1M2


Maschensatz fur M1: U0 = UR + U = R · I0 + U

Knotensatz fur K: I0 = I + C dUdt

Maschensatz fur M2: LdIdt = U

4! Achtung: Man beachte, dass in den Gleichungen die Große I0 auftritt, die keinem Energiespeicherzugeordnet ist. Daher muss sie aus den Gleichungen eliminiert werden, indem wir sie durch die Gleichungdes Knotensatzes ersetzen. Zum Schluss erhalten wir zwei DG fur die beiden Energiespeicher:

2.4 Einfache Beispiele mit L, C, R 21

⇒ U0 = R(I + C dUdt ) + U und LdI

dt = Ubzw.

⇒dUdt = ((U0 − U)/R− I)/C

dIdt = U/L

Algorithmus: U0 = A · sin(2πft)I = 0 und U = 0for t from tmin by dt to tmax

doU = U + ((U0 − U)/R− I)/C · dtI = I + U/L · dt...

end do

Typische Anwendung dieser Schaltung: Bezuglich der Spannung U wirkt die Schaltung als Bandpass:Hohe Frequenzen werden wegen der Wirkung von C, tiefe Frequenzen werden wegen der Wirkung vonL nicht durchgelassen. Nur ein Frequenzband um die mittlere Frequenz ω = 1√

LCkommt durch.

2.4.5 Parallelschwingkreis mit realen Spulen

C

R

L

I

U

K

U0

M1 M2

RL

RA


Maschensatz fur M1: U0 = UR + U = R · I0 + U

Knotensatz fur K: I0 = I + C dUdt + U

RA

Maschensatz fur M2: LdIdt + RLI = U

4! Achtung: Man beachte, dass wie im Falle des Parallelschwingkreises in den Gleichungen die GroßeI0 auftritt, die keinem Energiespeicher zugeordnet ist. Daher muss sie aus den Gleichungen eliminiertwerden, indem wir sie wieder durch die Gleichung des Knotensatzes ersetzen. Zum Schluss erhalten wirzwei DG fur die beiden Energiespeicher, die analog dem idealen Schwingkreis gelost werden.

⇒dUdt = ((U0 − U)/R− I − U/RA)/CdIdt = (U −RLI)/L


2.5 Aufstellen der DG fur Filterschaltungen

2.5.1 Aufstellen der DG fur einen Tiefpass (TP2PiCLC)

Gegeben ist der Tiefpass, der sich aus zwei Π-Gliedern zusammensetzt, TP2PiCLC:

Abb. 2.2. Tiefpass

Der Tiefpass besteht aus 5 Energiespeicher C1, C2, C3, L1, L2; diesen Energiespeichern werden 5 Zu-standsvariablen U1, U2, U3, I1, I2 zugeordnet. Systematisches Anwenden der Maschenregel und Kno-tenregel liefert die folgenden Gleichungen, wenn die physikalischen Gesetzmaßigkeiten UΩ = R · I(Ohmsches Gesetz), IC = C · U (Strom am Kondensator), UL = L · I (Spannungsabfall an derSpule) berucksichtigt werden:

M1: U0 = RI0 + U1

K1: I0 = C1U1 + I1 4! Achtung: I0 muss ersetzt werden, da es

M2: U1 = L1I1 + U2 keinem Energiespeicher zugeordnet ist!

K2: I1 = C2U2 + I2

M3: U2 = L2I2 + U3

K3: I2 = C3U3 + U3/R

Dies sind zunachst 6 Gleichungen, wobei Gleichung M1 keine DG darstellt. In dem zu losenden Systemdurfen als Variablen nur die 5 Zustandsvariablen U1, U2, U3; I1, I2 vorkommen, sonst keine. Also mussaus den Gleichungen M1 und K1 die Variable I0 eliminiert werden, da sie keinem Energiespeicher zuge-ordnet ist. Daher setzen wir Gleichung K1 in Gleichung M1 ein. Somit erhalten wir 5 DG 1. Ordnung:U1 = ((U0 − U1)/R− I1)/C1

I1 = (U1 − U2)/L1

U2 = (I1 − I2)/C2

I2 = (U2 − U3)/L2

U3 = (I2 − U3/R)/C3

Man beachte, dass pro DG nur eine Ableitung vorkommt und wir somit durch das Euler-Verfahren denfolgenden Algorithmus zur Losung des Systems erhalten

2.5 Aufstellen der DG fur Filterschaltungen 23

Algorithmus:U0 = . . . vorgegeben (s.u.)U1 = 0, I1 = 0, ... (Alle unbekannte Großen initialisieren)for t from tmin by dt to tmax

doU1 = U1 + ((U0 − U1)/R− I1)/C1 · dtI1 = I1 + (U1 − U2)/L1 · dtU2 = U2 + (I1 − I2)/C2 · dtI2 = I2 + (U2 − U3)/L2 · dtU3 = U3 + (I2 − U3/R)/C3 · dt...

end do

Bemerkungen:

(1) Man beachte, dass die DG in der Reihenfolge ihres Auftretens bei der physikalischen Modellierunggelost werden sollten, also von der Eingangsspannung zur Ausgangsspannung.

(2) Indem die Variablen im Algorithmus nicht mit Uneu1 und Ualt

1 bezeichnet werden, erspart man sichdie Umbenennung dieser Variablen und in die folgenden Gleichungen werden immer die aktuellen(also neu berechneten) Daten berucksichtigt.

Wir wahlen L1 = L2 = 1; C1 = C3 = 1; C2 = 2; R = 0.8 (kf = 1.314) (siehe Kapitel uber dieUbertragungsfunktion). Fur die numerische Rechnung wahlen wir ∆t = 0.005, tmax = 90.

(1.) Fur die Eingangsspannung setzen wir

U0(t) = 1. · sin(ωt) mit ω = 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0

und erhalten fur die Amplitude der Ausgangsspannung

Abb. 2.3. Ausgangssignal fur das Eingangssignal u0 = 1 · sin(ωt) mit ω = 1

ω 0.5 0.75 1.0 1.25 1.5 1.75 2.0UA 0.5 0.5 0.5 0.42 0.17 0.065 0.029

Man erkennt, dass die Ausgangsamplitude bei Frequenzen zwischen 1.25 und 1.5 drastisch abfallt.


(2.) Um einen Einschaltvorgang zu simulieren, setzen wir als Ein-

Abb. 2.4. Sprung S(t)

gangsspannung U0(t) = S(t). Dabei ist S(t) die Sprungfunktion.Fur negative Zeiten hat die Funktion den Wert 0 und springt beit = 0 auf den Wert 1. In Maple wird diese Funktion mit Heaviside(t)bezeichnet. Der Wert an der Stelle t = 0 muss explizit angegebenwerden.

Das zur Sprungfunktion gehorende Antwortsignal ist unten abgebildet.Da dieses die Reaktion des Systems auf die Sprungfunktion ist, nenntman das Antwortsignal die Sprungantwort:

Abb. 2.5. Sprungantwort

(3.) Wahlen wir als weiteres Eingangssignal (Eingangsspannung)

t

1/T

0 T

Abb. 2.6. Impulsstoß

einen Impulsstoß mit Breite T und Hohe 1/T . Damit regen wirdas System impulsartig an. Der Impulsstoß lasst sich in Maple mit1T (S(t)− S(t− T )) bzw.> U0(t) := 1/T*(Heaviside(t)-Heaviside(t-T));> Heaviside(0):=0;

realisieren. Die Ausgangsspannung nennt man zugehorig die Impulsantwort.

Abb. 2.7. Impulsantwort fur T = 0.5

2.5 Aufstellen der DG fur Filterschaltungen 25


Wir stellen die Differenzialgleichungen fur (1) einen Hochpass, (2) einen Tiefpass, (3) einen Bandpassund (4) eine Bandsperre auf und losen diese mit Maple. Die zugehorigen Maple-Worksheets konnendurch Anklicken gestartet und anschließend ausgefuhrt werden.

HP2TLCL Hochpass, der aus zwei T -Gliedern besteht.TP2PiCLC Tiefpass, der aus zwei Π-Gliedern besteht.BP2PiLCp Bandpass, der aus zwei Π-Gliedern besteht.BS2TLCp Bandsperre, die aus zwei T -Gliedern besteht.


2.6 Losen der DG fur den Tiefpass TP2PiCLC mit Maple

1. Aufstellen der DG> eq1:= R*(C1*dU1 + I1) = Ue-U1;

> eq2:= U1 = L1*dI1 + U2;

> eq3:= I1 = C2*dU2 +I2;

> eq4:= U2 = L2* dI2 + U3;

> eq5:= I2 = C3* dU3 +U3/R;

eq1 := R (C1 dU1 + I1 ) = Ue −U1

eq2 := U1 = L1 dI1 + U2

eq3 := I1 = C2 dU2 + I2

eq4 := U2 = L2 dI2 + U3

eq5 := I2 = C3 dU3 +U3R

2. Auflosen der DG nach den Ableitungen> sol:=solve(eq1,eq2,eq3,eq4,eq5, dU1,dU2,dU3,dI1,dI2);> assign(sol);

sol := dU2 = −−I1 + I2C2

, dI1 = −−U1 + U2L1

, dI2 = −−U2 + U3L2

, dU3 = −−I2 R + U3C3 R

,

dU1 = −R I1 −Ue + U1R C1

Parameter und Anfangsbedingungen> R:=0.8: C1:=1: C2:=C1*2: C3:=C1: L1:=1:L2:=L1:

> U1:=0: U2:=0: U3:=0: I1:=0: I2:=0:

Eingangsfunktion> Ue:=sin(w*t): w:=1.: #Wechselspannung

> Heaviside(0):= 0:

> Ue:=Heaviside(t): #Sprungfunktion

> T:=0.5: Ue:=Heaviside(t)-Heaviside(t-T):

3. Losen der DG mit dem Euler-Verfahren> tmax:=90:

> dt:=tmax/500:

> i:=0:

> for t from 0 by dt to tmax

> do i:=i+1:

> U1:=U1 + dt*dU1:

> I1 := I1 + dt* dI1:

> U2:= U2 + dt * dU2:

2.6 Losen der DG fur den Tiefpass TP2PiCLC mit Maple 27

> I2:= I2 + dt * dI2:

> U3:= U3 + dt * dU3:

> data1[i]:=[t,U3]:

> data2[i]:=[t,Ue]:

> od:

4. Darstellen der Losung> plot([seq(data1[n],n=1..i)], thickness=2);

Systemreaktion fur eine Wechselspannung als Eingangssignal bei w=1:

–0.4

–0.2

0

0.2

0.4

20 40 60 80

Systemreaktion fur eine Sprungfunktion als Eingangssignal:

0

0.1

0.2

0.3

0.4

0.5

20 40 60 80

Systemreaktion fur eine Impulsfunktion als Eingangssignal:

–0.02

0

0.02

0.04

0.06

0.08

20 40 60 80

3. Simulation von mechanischen Systemen mitMassenpunkten

Die Modellierung von mechanischen Systemen mit Massenpunkten gestaltet sich einfacher als die derelektrischen Filterschaltungen, da sie sich allgemein durch das Newtonsche Bewegungsgesetz beschrei-ben lassen:

Beschleunigungskraft = Masse * Beschleunigung= Summe aller angreifenden Krafte

m · −→r ′′ =∑n

i=1

−→Fi ,

wobei die Bewegungsgleichungen (=DG 2. Ordnung) komponentenweise angegeben werden. Z.B. furdie x-Komponente gilt

m · x′′ =n∑

i=1

(Fi)x = Fx

Diese DG 2. Ordnung transformiert man durch Einfuhrung der Geschwindigkeit (vx = x′) in ein Systemvon zwei DG 1. Ordnung:

v′x = ax =1m

Fx

x′ = vx

Wendet man auf jede einzelne DG das Euler-Verfahren an, so erhalt man den Algorithmus∣∣∣∣∣∣∣∣

vx = vx +1m

Fx dt

x = x + vx dt

t = t + dt

Im allgemeinen 3-dimensionalen Fall erhalt man somit 6 DG 1. Ordnung fur die 6 unbekannten Großenx(t), y(t), z(t), vx(t), vy(t), vz(t): ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣

vx = vx +1m

Fxdt

vy = vy +1m

Fydt

vz = vz +1m

Fzdt

x = x + vxdt

y = y + vydt

z = z + vzdt

t = t + dt

30 3 Simulation von mechanischen Systemen mit Massenpunkten

Beispiel: Schwingung eines 2-Federn-Masse-Systems

D1 D2

x(t)

m1

Abb. 3.1. 2-Federn-Masse-System

Eine Masse wird auf vertikaler Unterlage mit 2 Federn eingespannt. Die Ruhelage wird bei x = 0festgelegt. Lenkt man den Massenpunkt aus der Ruhelage nach rechts um x aus, so wirken auf dieMasse m sowohl die Ruckstellkraft der Feder D1 als auch von D2. Beide Krafte zeigen nach links, alsoin negative x-Richtung. Damit gilt nach dem Newtonschen Bewegungsgesetz

m · x′′ = −D1x−D2x

Berucksichtigt man noch eine Reibungskraft, die der Geschwindigkeit entgegen wirkt, FR ∼ −v = −x′,erhalt man

x′′(t) = −D1 + D2

m· x(t)− γ · x′(t).

Dies fuhrt auf den folgenden Algorithmus zum Losen des DG Systems 1. Ordnung

v = v + (−(D1 + D2)/m · x− γ · v) dt

x = x + v dt

t = t + dt

mit den Anfangsbedingungen x(0) = x0 und v(0) = 0. Dieser Algorithmus zusammen mit der numeri-schen Losung ist im ersten Teil des zugehorigen Maple-Worksheet umgesetzt.

3 Simulation von mechanischen Systemen mit Massenpunkten 31

Beispiel: Schwingung eines 2Massen-3Federn-Systems

D1 D2 D3

x (t)2x (t)1

m2

m1

Abb. 3.2. 2Massen-3Federn-Systems

Zwei Massen m1 und m2 werden gemaß obigem Bild durch 3 Federn (mit Federkonstanten D1, D2, D3)gekoppelt. Die Auslenkungen aus der Ruhelage werden mit x1(t) fur Masse m1 und x2(t) fur Masse m2

bezeichnet. Lenkt man die Massen m1 und m2 aus, so wirken die Federn D1 und D3 auf die absolutenAuslenkungen x1 bzw x2. Fur die Feder D2 geht jedoch die relative Auslenkung x1 − x2 ein. Es geltendie Bewegungsgleichungen

m1 · x′′1(t) = −D1 · x1(t)−D2(x1(t)− x2(t))− γ · x′1(t)m2 · x′′2(t) = −D3 · x2(t)−D2(x2(t)− x1(t))− γ · x′2(t)

wenn man eine lineare geschwindigkeitsabhangige Reibungskraft berucksichtigt. Zur numerischen Lo-sung fuhrt man die beiden Geschwindigkeiten v1(t) = x′1(t) und v2(t) = x′2(t) ein und erhalt ein DGSystem 1. Ordnung der Form

v′1(t) = −D1 + D2

m1x1(t) +

D2

m1x2(t)−

γ

m1v1(t)

x′1(t) = v1(t)

v′2(t) = −D3 + D2

m2x2(t) +

D2

m2x1(t)−

γ

m2v2(t)

x′2(t) = v2(t).

Dieses DG System fuhrt auf den folgenden Algorithmus

x1=0.1, x2=0, v1=0, v2=0 (Initialisierung)

for t from 0 by dt to tmax dov1 = v1 + dt*(-(D1+D2)/m1*x1 + D2/m1*x2 - gamma/m1*v1)v2 = v2 + dt*(-(D3+D2)/m2*x2 + D2/m2*x1 - gamma/m2*v2)

x1 = x1 + dt * v1x2 = x2 + dt * v2

end do

32 3 Simulation von mechanischen Systemen mit Massenpunkten

Auslenkung Masse 1

–0.1

0.05

0

0.05

0.1

20 40 60 80 100

Auslenkung Masse 2

–0.1

0.05

0

0.05

0.1

20 40 60 80 100

Man erkennt als Schwingungsform eine Schwebung: Masse m1 schwingt maximal; Masse m2 ist zu-nachst noch in Ruhe. Uber die Feder D2 wird aber Energie von Masse m1 auf die Masse m2 ubertragen,diese beginnt zu schwingen. In dem Maße wie die Schwingung von Masse m2 zunimmt, nimmt gleich-zeitig die Energie von m1 ab bis m1 sogar zur Ruhe kommt und m2 maximal schwingt. Dann kehrtsich der Prozess um.

Dieser Algorithmus zusammen mit der numerischen Losung ist im zweiten Teil des zugehorigen Maple-Worksheet umgesetzt.

4. Simulation geladener Teilchen inelektromagnetischen Feldern

Die Bewegung von elektrisch geladenen Teilchen (Masse m, Ladung q) in elektrischen und magnetischenFeldern ist bestimmt durch die Lorentz-Kraft

−→F L = q(

−→E +−→v ×

−→B ).

Zur Bestimmung der magnetischen Kraft berechnen wir

−→v ×−→B =

∣∣∣∣∣∣ex vx Bx

ey vy By

ez vz Bz

∣∣∣∣∣∣ =

vyBz − vzBy

−vxBz + vzBx

vxBy − vyBx

und erhalten als Bewegungsgleichungen

m

x(t)y(t)z(t)

′′

= q(

Ex

Ey

Ez

+

vyBz − vzBy

−vxBz + vzBx

vxBy − vyBx

).

Spezialfall: Das Magnetfeld besitzt nur eine z-Komponente (Bx =

V0

B=B eZ0

Abb. 4.1. Magnetfeld

0, By = 0, Bz = B0), d.h. es wirkt nur in z-Richtung. Dann erhaltman

m · x′′(t) = q · Ex + q ·B0 · vy

m · y′′(t) = q · Ey − q ·B0 · vx

m · z′′(t) = q · Ez .

Fuhrt man wie im mechanischen Fall die Geschwindigkeit als zusatzliche Große ein (vx, vy, vz) gilt

v′x =q

mEx +

q

m·B0 · vy

v′y =q

mEy −

q

m·B0 · vx

v′z =q

mEz

x′ = vx

y′ = vy

z′ = vz

34 4 Simulation geladener Teilchen in elektromagnetischen Feldern

Das fuhrt auf den Algorithmus

vx = vx + (q

mEx +

q

m·B0 · vy) dt

vy = vy + (q

mEy −

q

m·B0 · vx) dt

vz = vz +q

mEz dt

x = x + vx dt

y = y + vy dt

z = z + vz dt

t = t + dt

1. Fall: Bewegung eines geladenen Teilchens nur im Magnetfeld. Dann ist das elektrische Feld kom-

plett Null,−→E = 0. Die Bewegung ist dann in der (x,y)-Ebene eine Kreisbahn. Die Geschwindigkeit in

z-Richtung bleibt konstant, da v′z = 0:

–1–0.5

00.5

1

00.5

11.5

2

00.5

11.5

22.5

3

(a) Flugbahn des Teilchens

0

0.5

1

1.5

2

5 10 15 20 25 30

(b) Nur die y-Komponente

2. Fall: Bewegung eines geladenen Teilchens in elektrischen und magnetischen Feldern. Die Kreisbahnwird dann durch eine Beschleunigung uberlagert.

020

4060

80100

120

160–120

–80–40

0

–40

–30

–20

–100

Abb. 4.2. Flugbahn des Teilchens

5. Signal- und Systemanalyse mit der FFT

5.1 Grundlegende Formeln

Motivation: Gegeben ist ein Zeitsignal, welches nur in Form von diskreten Werten (t0, f0), (t1, f1),. . . , (tn−1, fn−1) vorliegt. Gesucht sind die Frequenzen und zugehorigen Amplituden, die in diesemSignal enthalten sind.

Abb. 5.1. Zeitsignal

Kontinuierliche Signale: Die Formeln fur die Fourieranalyse (Fourierreihen fur periodische Signale,Fouriertransformation fur beliebige Zeitsignale) setzen voraus, dass die Funktion f(t) als kontinuierlichesSignal in Form einer Funktionsvorschrift vorliegt. Gehen wir zur Wiederholung der Formeln zunachstvon einem T-periodischen Signal aus.

Abb. 5.2. T-periodisches Signal

Dann lassen sich die komplexen Fourierkoeffizienten cm der Fourierreihe bestimmen durch

cm =1T

∫ T

0

f(t)e−imω0tdt; mit w0 =2π

T.

36 5 Signal- und Systemanalyse mit der FFT

cm sind die komplexen Fourierkoeffizienten; sie geben bis auf den Faktor 2 an mit welcher komplexenAmplitude (Betrag und Phase) die Frequenz ωm = mω0 im Signal enthalten ist. In einem periodischenSignal sind also nur die diskreten Frequenzen

ω0, 2ω0, 3ω0, . . . ,mω0

enthalten, der Frequenzabstand betragt

∆ω = ω0 =2π

T.

Diskretes Signal: Kommen wir nun zu dem Fall, dass

Abb. 5.3. diskretes Signal

f in Form einer Abtastung vorliegt. Dann gibt es zuf nur diskrete Funktionswerte zu den gegebenen Zei-ten t0, t1, . . . , tN−1, aber keine explizite Funktions-vorschrift. Wir setzen f T-periodisch fort und kon-nen die komplexen Fourierkoeffizienten naherungswei-se bestimmen durch

cm =1T

∫ T

0

f(t)e−imω0tdt ≈ 1T

N−1∑n=0

f(tn)e−imω0tn∆t.

Mit tn = n∆t = n TN und ω0 = 2π

T folgt weiter

cm ≈ 1T

N−1∑n=0

fne−im 2πT n T

NT

N=

1N

N−1∑n=0

fne−imn 2πN︸︷︷︸

DFT

Man nennt ∑N−1n=0 fne−imn 2π

N

die diskrete Fouriertransformation (DFT). Sie gibt bis auf einen Normierungsfaktor ( 2N ) )1 an, mit

welcher Amplitude die Frequenz mω0 im Signal naherungsweise enthalten ist. Wieder erhalt man eindiskretes Spektrum mit den Frequenzen ω0, 2ω0, 3ω0, . . . mit dem Frequenzabstand ∆ω = 2π

T .

Zugang uber die Fouriertransformation: Wir vergleichen die eben erhaltene Formel fur die DFT mitder Fouriertransformation:

F (ω) =∫ ∞

−∞f(t)e−iωtdt =

∫ T

0

f(t)e−iωtdt

≈N−1∑n=0

f(tn)e−iωtn∆t

Dies ist zunachst eine kontinuierliche Funktion in ω. Speziell fur die diskreten ω-Werte ωm = mω0 giltaber

1 Man beachte, dass die komplexen Fourierkoeffizienten bis auf den Faktor 2 dem Amplitudenspektrumder periodischen Funktion entsprechen (siehe Mathematik III) Am =

√a2

m + b2m = 2|cm|

5.2 Eigenschaften der DFT 37

F (mω0) ≈N−1∑n=0

fne−im 2πT n T

NT

N

≈ T

N

N−1∑n=0

fne−imn 2πN︸︷︷︸

DFT

F (mω0) stimmt also bis auf den Faktor TN mit der DFT uberein:

F (mω0) = Tcm.

Im Folgenden kurzen wir die Formel fur die DFT mit dem Symbol Fm ab. Diese zwar etwas laxe abereindeutige Schreibweise soll daran erinnern, dass bis auf den Normierungsfaktor 2

N der Betrag |Fm| dem

reellen Amplitudenspektrum der Funktion f entspricht; Fm aber auch bis auf einen NormierungsfaktorTN dem Wert der Fouriertransformation an der Stelle mω0 entspricht:

Fm :=N−1∑n=0

fn e−imn 2πN (DFT) (5.1)

Fm wird im Folgenden auch als das Spektrum bezeichnet.

5.2 Eigenschaften der DFT

(1) Fur m ≥ N wiederholen sich die Spektren:

Fm+N =N−1∑n=0

fne−i(m+N)n 2πN =

N−1∑n=0

fne−imn 2πN · e−in2π︸︷︷︸

=1

= Fm

D.h. fur m ≥ N erhalt man keine neue Frequenzinformation, sondern das Spektrum wiederholt sichperiodisch.


(2)Ist f ein reelles Signal, dann liegt das Spektrum symmetrisch (bis auf komplexe Konjugation)zu N

2 :

FN−k =N−1∑n=0

fne−i(N−k)n 2πN =

N−1∑n=0

fneikn 2πN · e−in2π︸︷︷︸

=1

=N−1∑n=0

fne−ikn wπN = Fk

D.h. fur m ≥ N2 erhalt man keine neue Frequenzinformation, sondern die Amplituden sind komplex

konjugiert symmetrisch um N2 .

(3)

Folgerung: Zu N Messwerten f0, f1, . . . , fN−1 kann man durch die DFT hochs-tens die Amplituden zu den Frequenzen ω0, . . . , ωN

2ablesen, wobei die hochste

noch erkennbare Frequenz ωN2

= N2 ω0 = N

22πT ist.

Die Abtastfrequenz betragt

ωa = 2πfa = 2π1

∆t= 2π

N

T= N

2π

T⇒ ωa = 2ωN

2

Diese Folgerung besagt, dass die Abtastfrequenz ωa mindestens doppelt so groß sein muss als diemaximal im Signal enthaltene Frequenz, um sie durch die Analyse zu erkennen (=Abtasttheorem vonShannon).

4! Achtung: Man beachte, dass die Abtastfrequenz ωa echt großer sein muss als 2ωmax, daωa = 2ωmax nicht ausreicht, wie das folgende Beispiel zeigt:

Beispiel: Wahlt man als Signal sin(2πt), so ist die Frequenz ω = 2π. Tastet man das Signal mit derdoppelten Frequenz ωa = 4π ab, so entspricht dies den Funktionswerten zu den Zeit mit ∆t = 1

2 .

sin(2πt) ⇒ ωmax = 2π ⇒ ωa = 4π ⇒ ∆t =2π

ωa=

12

Damit erhalt man als abgetastete Funktionswerte nur die Werte Null, welche das Signal naturlich nichtreprasentieren.

Abb. 5.4. Falsche Abtastung eines Sinussignals

5.3 Berechnung der DFT 39

(4) Die inverse DFT erhalt man uber die Formel

f(tn) =N−1∑k=0

ckeikn 2πN =

1N

N−1∑k=0

Fk eikn 2πN .

Man beachte, dass sich die inverse DFT im Wesentlichen mit der gleichen Formel berechnet wie dieDFT selbst, denn es gilt

fn =1N

N−1∑k=0

Fke−ikn 2πN

D.h. man bildet von Fk die DFT und nimmt davon das komplex konjugierte.

Da Fk komplexe Großen darstellen, programmiert man die DFT gleich allgemein fur komplexe Abtast-werte, obwohl man sie fur die Analyse reeller Signale so allgemein nicht benotigte. Bei der Verwendungvon Programmroutinen zur Berechnung der DFT muss man daher beachten, dass im Falle von reellenSignalen der Imaginarteil auf Null gesetzt werden muss!

5.3 Berechnung der DFT

(1.) Einfacher Algorithmus zur Berechnung der DFT: Gegeben sind N diskrete Messwerte

f0, f1, . . . , fN−1. Gesucht ist Fm, d.h. die komplexe Amplitude zur Frequenz ωm:

Fm =N−1∑n=0

fne−imn 2πN =

N−1∑n=0

fn (cos(mn2π

N)− i · sin(mn

2π

N))

=N−1∑n=0

fn cos(mn2π

N)− i

N−1∑n=0

fn sin(mn2π

N).

Zur Berechnung der Amplitude Fm fur ωm wird folgende Schleife abgearbeitet:re := 0;im := 0;for n from 0 to N-1do

re := re + f[n] * cos( m*n*2*Pi/N ):im := im + f[n] * sin( m*n*2*Pi/N ):

end do:Am := 2 * sqrt(re^2 + im^2):

Fur jede Frequenz muss obige Schleife mit anderem m durchlaufen werden. Bei 1000 Messwerten mus-sen fur die Bestimmung einer Amplitude 1000 Sinus- und Kosinusausdrucke ausgewertet werden →sehr, sehr zeitaufwendig.


(2.) Schneller Algorithmus zur Berechnung der DFT: FFTDie FFT ist ein schneller Algorithmus, um

Fm =N−1∑n=0

fne−imn 2πN fur m = 0, . . . , N − 1

auszurechnen. Voraussetzung hierfur ist, dass

N = 2k

Zusammenfassung

(1) Sind f0, f1, . . . , fN−1 gegebene Messwerte mit N = 2k, dann erhalt man mit der FFT dieAmplituden zu den Frequenzen

ω0, ω1, . . . , ωN2,

d.h. bis zur maximalen Frequenz ωN2

= N2 ω0 = N

22πT . Der Frequenzabstand betragt

∆ω = 2πT .

(2) Um die Frequenzauflosung zu erhohen (d.h. kleineres ∆ω), vergroßert man T .

(3) Um die maximale Frequenz zu vergroßern wird bei festem T die Anzahl der Abtastungen Nerhoht.

(4) Die DFT bzw. FFT liefern komplexwertige Großen. Um das reelle Amplitudenspektrum zuerhalten, geht man zum Betrag uber.

Am = 2√

Re(Fm)2 + Im(Fm)2

(5) Die FFT funktioniert nur, wenn die Anzahl der Messwerte eine Potenz von 2 darstellt: N = 2k.

(6) Um vorhandene Programmroutinen zur FFT auf reellwertige Signale anzuwenden, muss derImaginarteil auf Null gesetzt werden.

5.4 FFT mit Maple 41

5.4 FFT mit Maple

Im folgenden Kapitel wird die Signal- und Systemanalyse mit Maple durchgefuhrt und die Ergebnissediskutiert. Es zeigt sich, dass erst durch die geeignete Wahl der Abtastparameter eine aussagekraftigeInterpretation moglich wird.

5.4.1 Signalanalyse

Im Worksheet FFT wird eine Signalanalyse mit Maple durchgefuhrt. Es werden dabei die folgendenAspekte diskutiert:

(1) FFT eines Sinussignals mit Maple: Wie hangt der Spektralbereich von der Anzahl der Abtastpunkteab?

(2) Analyse mit dem Hanning-Filter.

(3) Die Prozedur plot fft mit der Option Hanning-Fenster.

(4) Maximal erfassbare Frequenz.

(5) Frequenzauflosung.

(6) Analyse von Messdaten.

5.4.2 Systemanalyse

Im Worksheet sys ana wird eine Frequenzanalyse eines Tiefpasses, der aus zwei Π-Gliedern zusammen-gesetzt ist (TP2PiCLC), durchgefuhrt.

6. Signal- und Systemanalyse mit Maple

6.1 FFT eines Sinussignals mit Maple: Wie hangt der Spektralbereichvon der Anzahl der Abtastpunkte ab?

Wie der Spektralbereich von der Anzahl der Abtastpunkte abhangt, sei am Beispiel der Funktionsin(2π t) demonstriert:> restart:

> f:=t -> evalf(sin(2*Pi*t)): #Definition der Funktion

> plot(f(t),t=0..10.1,color=black, thickness=2);

–1

–0.5

0

0.5

1

2 4 6 8 10t

(1.) Abtastung der Funktion: Die Abtastung der Funktion f(t) erfolgt im Intervall [0..10.1] an denStellen t=i*dt mit dt=T/N.> m:=8: N:=2m; #Festlegung der Abtastrate

> T:=10.1: dt:=evalf(T/N): #Festlegung des Zeitintervalls

> for i from 1 to N

> do fd[i] := f((i-1)*dt):

> imd[i]:= 0:

> end do:

N := 256

44 6 Signal- und Systemanalyse mit Maple

(2.) Berechnung der DFT mit der FFT: Mit dem Maple-Befehl FFT wird das komplexe Spektrum derFunktion berechnet> FFT(m, fd, imd);

256

(3.) Graphische Darstellung des Spektrums: Die Aufbereitung der Daten erfolgt dadurch, dass wirvon den komplexen Großen (fd, imd) zum Betrag ubergehen und sie als Funktion der Frequenzenw=0..N/2*w0 darstellen. Der Skalierungsfaktor, mit dem die Transformierten multipliziert werden, ist2/N.> for i from 1 to N/2

> do daten[i]:= [2*3.14/T*(i-1), 2/N*sqrt(fd[i] 2+imd[i] 2)]:

> end do:

Fur den plot-Befehl stellen wir die Daten als Liste von Wertepaaren [[x1, y1], [x2, y2], ..., [xn, yn]] bereit:> plot([seq(daten[i],i=1..N/2)], labels=[”w”,” ”],color=’red’, thickness=2);

0

0.2

0.4

0.6

0.8

20 40 60 80w

Aufgabe: Man variiere die Abtastrate N=64, 256, 1024, 2048 und vergleiche die Ergebnisse. Anschlie-ßend stelle man fur die Falle N=256, 1024, 2048 nur den Bereich w=0..20 dar.

Zusammenfassung: Durch Variation der Abtastrate N andert man bei festem T nur w max; dieFrequenzauflosung bleibt dabei erhalten; insbesondere wird die Linie bei w=6.28 durch Vergro-ßerung von N nicht scharfer dargestellt!

Aufgabe: Man variiere bei festem N=256 die Abtastzeit T : T=10.1, T=20.1, T=30.1. Wie andertsich die Breite der Spektrallinie? Man wahle T=10 und N=256. Warum erhalt man mit T=10 genaueine Linie mit der DFT? → Hanning-Fenster → plot fft

6.2 Das Hanning-Fenster 45

6.2 Das Hanning-Fenster

Problem: Wenn man periodische Signale in einem Zeitfenster betrachtet, welches nicht mit der Periodezusammenfallt, dann schneidet man das Signal mit einem falschen T ab. Denn wahlt man z.B. fur dieSinusfunktion nicht ein Vielfaches von 2 Pi> T:=8.4:

> plot(sin(t), t=0..T);

–1

–0.5

0

0.5

1

2 4 6 8t

so wird das Signal fur die Analyse nun periodisch mit der Periode T=8.4 fortgesetzt. Damit kommt esaber genau an der Stelle t=T zu einer Sprungstelle im periodisch fortgesetzten Signal.> signal:=piecewise(0<t and t<T, sin(t), T<t and t<2*T, sin(t-T), 0):

> plot(signal,t=0..2*T, discont=true);

–1

–0.5

0

0.5

1

2 4 6 8 10 12 14 16t

Die Fourierkoeffizienten einer periodischen Funktion mit Sprungstelle verhalten sich wie ∼1/n, d.h.Oberfrequenzen dominieren in dem Signal, obwohl diese Frequenzen in dem ursprunglichen Sinus nichtvorhanden sind. Außerdem besitzt das Signal im Intervall [0, T] einen Gleichanteil, da der lineareMittelwert ungleich Null ist> 1/T*int(sin(t), t=0..T);

.1808676968

Ausweg: Um beide kunstlichen Effekte (=Abschneide-Effekte) zu eliminieren, multipliziert man das

Signal mit einer Funktion, die bei t=0 und t=T Null ist, namlich mit der Funktion 1− cos( 2 π tT )


> plot([sin(t), 1-cos(2*Pi/T*t)], t=0..T);

–1

–0.5

0

0.5

1

1.5

2

2 4 6 8t

Im unten angegebenen Bild ist sowohl das Signal (rot), die Funktion 1 − cos( 2 π tT ) (grun) sowie das

amplitudenmodulierte Signal (blau) eingezeichnet. Die FFT wird nun anstatt von dem ursprunglichenSignal nun von dem amplitudenmodulierten Signal berechnet.> plot([sin(t),1-cos(2*Pi/T*t),sin(t)*(1-cos(2*Pi/T*t))], t=0..8.2, color=[red, green,blue], thickness=[1,1,3]);

–2

–1

0

1

2

2 4 6 8t

Man bezeichnet dieses Vorgehen als Hanning-Filter (Siehe auch den nachsten Abschnitt: Die Prozedurplot fft und Hanning-Fenster). Die Modulationsfunktion ist so gewahlt, dass sie das ursprungliche Signalan den Enden glatt auf Null druckt und der Mittelwert der Funktion 1 ergibt.> 1/T*int((1-cos(2*Pi/T*t)),t=0..T);

1.000000000

6.3 Die Prozedur plot fft mit der Option Hanning 47

6.3 Die Prozedur plot fft mit der Option Hanning

Wir tasten das Signal sin(2π t) nochmals ab und verwenden die Prozedur plot fft, um eine Frequenz-analyse durchzufuhren> f:=t -> evalf(sin(2*Pi*t)): #Definition der Funktion

> plot(f(t),t=0..10.5,color=black, thickness=2);

–1

–0.5

0

0.5

1

2 4 6 8 10t

(1.) Abtastung der Funktion: Die Abtastung der Funktion f(t) erfolgt im Intervall [0..10.5] an denStellen t = i ∗ dt mit dt = T/N .> m:=8: N:=2m; #Festlegung der Abtastrate

> T:=10.5: dt:=evalf(T/N): #Festlegung des Zeitintervalls

> for i from 1 to N

> do fd[i] := f((i-1)*dt):

> end do:

(2.) Berechnung der DFT und die graphische Darstellung des reellen Spektrums mit der Prozedurplot fft: Die Prozedur plot fft wird aus dem lokalen Verzeichnis eingelesen. Dazu kopiert man sich dieDatei prozeduren.m in das lokale Verzeichnis und liest die Prozedur ein mit> read(prozeduren.m);

Ist die Datei in einem anderen Verzeichnis als das Maple-Worksheet, das man gerade verwendet, dannmuss man den vollen Pfad angeben mit> read(E://verz1//verz2//prozeduren.m);

Die Transformation und die Aufbereitung der Daten erfolgt nun durch > plot fft(fd, m, T);

1. Argument: Abgetastete Werte als Liste [f1, f2, ..., fN]

2. Argument: m: 2ˆm ist die Anzahl der Abtastungen

3. Argument: Abtastzeit

Das berechnete Spektrum ist richtig zu interpretieren, indem eine dominierende Frequenz, namlichw=6.2, im Signal enthalten ist. Das berechnete Spektrum ist aber falsch zu interpretieren, indem zum


0

0.1

0.2

0.3

0.4

0.5

0.6

A

10 20 30 40 50 60 70w

Abb. 6.1. FFT ohne Hanning-Filter

einen das Signal aber auch einen Gleichstromanteil besitzt, da die Amplitude bei w=0 nicht Null ergibt.Und zum anderen auch hohere Frequenzen im Signal vorkommen.

Beides ruhrt von dem falschen Abschneiden der periodischen Funktion her. Durch das falsche Abschnei-den kommt es bei der Stelle T dann zu einem Sprung → hohe Frequenzen kommen in der Fourierreihevor (siehe Mathe III) und das Signal hat auch einen Offset → es gibt einen Gleichspannungsanteil →Amplitude bei w=0. Legt man uber das Zeitfenster nun ein Hanning-Fenster (die Funktion bugelt dieAbtastwerte bei t=0 und bei t=T zu Null) treten beide oben beschriebenen Effekte nicht mehr auf.Das Hanning-Fenster ist als Option in der plot fft Routine berucksichtigt:

> plot fft(fd, m, T, Hanning);

0

0.2

0.4

0.6

0.8

A

10 20 30 40 50 60 70w

Abb. 6.2. FFT ohne Hanning-Filter

Diskussion: Durch Verwendung des Hanning-Fensters werden Frequenzen, die kunstlich in derAnalyse durch das Abschneiden des Signals bzw. dem endlichen Zeitfensters erscheinen, eliminiert.

6.4 Maximal erfassbare Frequenz 49

6.4 Maximal erfassbare Frequenz

Gegeben ist ein Signal, das sich aus mehreren Frequenzen zusammensetzt. Gesucht wird die Abtastung,bei dem die großte im Signal enthaltene Frequenz noch aufgelost wird. Vorgegeben ist eine Uberlagerungvon Sinusfunktionen mit den Frequenzen 1 Hz, 2 Hz, 3 Hz und 4 Hz.> f:=t -> evalf(sin(2*Pi*1*t) + sin(2*Pi*2*t) + sin(2*Pi*3*t) + sin(2*Pi*4*t) ):

> plot(f(t),t=0..10,color=black, thickness=2);

–3

–2

–1

0

1

2

3

2 4 6 8 10t

(1.) Abtastung der Funktion f(t) an den Stellen t=i*dt mit dt=T/N.> m:=6: N:=2 m:

> T:=10.1: dt:=evalf(T/N):

> for i from 1 to N

> do fd[i] := f((i-1)*dt):

> end do:

(2.) Berechnung der FFT> read(prozeduren.m):

> plot fft(fd, m, T);

0.2

0.4

0.6

0.8

1

A

0 5 10 15w


Diskussion: Aus dem Diagramm entnimmt man zwar, dass im Signal 4 Frequenzen vorkommen, aberwie man aus der Definition des Signals erkennen kann ist eine davon falsch dargestellt. Letztendlich liegtes daran, dass die hochste im Signal auftretende Frequenz nicht richtig abgetastet wird und daher in denkleineren Frequenzbereich gespiegelt wird ( Aliasing-Effekt). Um dies zu verdeutlichen/zu verstehen,fuhre man die folgenden Aufgaben durch:

Aufgabe: Man variiere die Abtastrate N=64, 128, 256, 512, 1024 und vergleiche die Ergebnisse mitN=256. → Aliasing-Effekt

Zusammenfassung: Durch Variation der Abtastrate N andert man bei festem T nur wmax.Ist N zu klein gewahlt, werden nicht mehr alle Frequenzen des Signals richtig erfasst. Nach demAbtasttheorem muss die Abtastfrequenz wa > 2ws = 2∗25 1/s gewahlt werden, um die Frequenz25 1/s noch zu erfassen. Aus wmax = 2π/TN folgt fur eine Auflosung der Frequenz von 25 1/s:N > 80.

6.5 Frequenzauflosung

Gegeben ist wieder ein Signal, das sich aus mehreren Frequenzen zusammensetzt. Gesucht ist derAbtastungsparameter, bei dem alle Frequenzen getrennt aufgelost werden. Diese Untersuchung sei amBeispiel der Uberlagerung von 4 Sinusfunktionen mit den Frequenzen 1.0 Hz, 1.2 Hz, 2.3 Hz und 2.7Hz durchgefuhrt.> f:=t->evalf(sin(2*Pi*1.0*t) + sin(2*Pi*1.2*t) + sin(2*Pi*2.3*t) + sin(2*Pi*2.7*t)):

> plot(f(t),t=0..15,color=black, thickness=2);

–3

–2

–1

0

1

2

3

2 4 6 8 10 12 14t

(1.) Abtastung der Funktion f(t) an den Stellen t=i*dt mit dt=T/n.> m:=8: N:=2 m;

> T:=5.1: dt:=evalf(T/N):

6.6 Analyse von Messdaten 51

> for i from 1 to N

> do fd[i] := f((i-1)*dt):

> end do:

(2.) Berechnung der FFT> read(prozeduren.m):


0

0.2

0.4

0.6

0.8

1

A

20 40 60 80 100 120 140w

Diskussion: Aus dem Diagramm entnimmt man, dass im Signal 3 Frequenzen vorkommen, aber wieman aus der Definition des Signals erkennen kann ist dies falsch, da 4 Frequenzen festgelegt werden.Letztendlich liegt es daran, dass zwei Frequenzen unter einem Peak liegen, d.h. die Auflosung nichthoch genug ist. Um dies zu verdeutlichen/zu verstehen, fuhre man die folgenden Aufgaben durch:

Aufgabe: Man variiere T und N so, dass eine gute Auflosung der 4 Frequenzen gewahrleistet ist.

Zusammenfassung: Um die Auflosung des Spektrums zu erhohen, d.h. einzelne Spektrallinieneines Signals besser auflosen zu konnen, muss die Messdauer T erhoht werden. Die Anzahl derAbtastungen N hat keinen Einfluss auf die Auflosung!

6.6 Analyse von Messdaten

Die plot fft-Routine kann auch verwendet werden, um Messdaten zu analysieren: In der Textdateidaten1.txt sind Messdaten abgespeichert. Die erste Spalte besteht aus den Zeitwerten ti, die zweiteSpalte aus den abgetasteten Funktionswerten. Diese Daten kann man z.B. auch mit einem Texteditoranschauen. Mit dem readdata-Befehl kann man diese Daten in Maple einlesen und einer Variablenzuweisen. Es ist beim Einlesen zu beachten, dass sich die Datei in dem Verzeichnis befindet, aus demman Maple gestartet hat. Dann kann man einen lokalen Namen verwendenliste:=readdata(daten1.txt, 2):


Ansonsten muss der globale Dateiname angegeben werden, z.B.liste:=readdata(E://verz1//verz2//daten1.txt, 2):In jedem Fall aber muss als zweites Argument eine 2 gesetzt werden. Die 2 steht fur eine 2-spaltigeDatei.> liste:= readdata(”daten1.txt”,2):

Mit> liste[10];

[0.001757812500, −2.580593604]

kann man z.B. die 10. Zeile ansteuern: die erste Zahl steht fur den Zeitpunkt, die zweite Zahl fur denabgetasteten Funktionswert zu diesem Zeitpunkt. Mit> N:=nops(liste);

N := 256

lasst man sich die Anzahl der Messdaten angeben. Die Messzeit kann man damit bestimmen, indemman in die N.-te Zeile geht (liste[N]) und von dieser letzten Zeile die erste Zahl> T:= liste[N][1];

T := 0.04980468750

(1.) Graphische Darstellung der Messdaten> plot(liste);

–10

–5

0

5

10

0.01 0.02 0.03 0.04 0.05

(2.) Anwenden der plot fft-Routine: Man beachte, dass bei der plot-fft-Routine ja nur die diskretenFunktionswerte benotigt werden und nicht die Paare (t i, f i). Daher bestimmt man zuerst fd als diezweite Spalte der Liste und setzt bei N=256 Daten m auf 8> fd:= [seq( liste[i][2],i=1..N )]:

> m:=8:

> read(prozeduren.m):


6.6 Analyse von Messdaten 53

0

2

4

6

8

A

5000 10000 15000w

Man kann aus dem Diagramm 3 Frequenzen ablesen: w1 = 5000, w2 = 10000, w3 = 12000.


6.7 System-Analyse mit Maple

Die DFT (und deren schnelle Umsetzung FFT) wird messtechnisch sowohl zur Analyse von Signalen (→Signalanalyse) als auch zur Bestimmung des Ubertragungsverhaltens von Systemen (→ Systemanalyse)intensiv verwendet. Die Signalanalyse mit der FFT stellt sich schematisch dar als

Signal —-> FFT —-> Spektrum

Die Systemanalyse mit der FFT lasst sich schematisch im folgenden Diagramm darstellen:

Signal —-> SYSTEM —-> Systemreaktion

↓ ↓

Spektrum des Spektrum desEinganssignals Ausgangssignals

Die Systemfunktion (=Ubertragungs- oder Transferfunktion) ist das Verhaltnis vom Spektrum des Aus-gangssignals zum Spektrum des Eingangssignals

H(w) =G(w)F(w)

=FFT(g)FFT(f)

wenn f(t) das Eingangssignal und g(t) das zugehorige Ausgangssignal bzw. G(w) das Spektrum desAusgangs- und F(w) das Spektrum des Eingangssignals ist.

Beispiel: Frequenzanalyse eines Tiefpasses. Gegeben ist ein Tiefpass, der aus zwei Pi-Gliedern zu-sammengesetzt ist (TP2PiCLC).

(1.) Bestimmung der diskreten Impulsantwort: Um die diskreten Werte der Impulsantwort zu er-halten, losen wir numerisch die das Netzwerk beschreibenden DG. Die numerischen Werte der Impul-santwort entsprechen den abgetasteten Werten der Ausgangsspannung U3 (=Ausgangssignal).

6.7 System-Analyse mit Maple 55

Aufstellen der DG> restart:

> eq1:= R*(C1*dU1 + I1) = Ue-U1;

> eq2:= U1 = L1*dI1 + U2;

> eq3:= I1 = C2*dU2 +I2;

> eq4:= U2 = L2* dI2 + U3;

> eq5:= I2 = C3* dU3 +U3/R;

eq1 := R (C1 dU1 + I1 ) = Ue −U1

eq2 := U1 = L1 dI1 + U2

eq3 := I1 = C2 dU2 + I2

eq4 := U2 = L2 dI2 + U3

eq5 := I2 = C3 dU3 +U3R

und Auflosen der DG nach den Ableitungen> sol:=solve(eq1,eq2,eq3,eq4,eq5, dU1,dU2,dU3,dI1,dI2);> assign(sol);

sol := dU1 = −R I1 −Ue + U1R C1

, dU3 =I2 R−U3

C3 R,

dI2 =U2 −U3

L2, dU2 =

I1 − I2C2

, dI1 =U1 −U2

L1

Die Parameter der Bauelemente und die Anfangsbedingungen werden festgelegt> R:=0.8: C1:=1: C2:=C1*2: C3:=C1: L1:=1:L2:=L1:

> U1:=0: U2:=0: U3:=0: I1:=0: I2:=0:

Fur die Eingangsfunktion wahlen wir eine Naherung an den Delta-Impuls> t0:=0.5:

> Heaviside(0):=1:

> Ue:=(Heaviside(t)-Heaviside(t-t0))/t0:

> plot(Ue, t=0..90, numpoints=500, thickness=2);

0

0.5

1

1.5

2

20 40 60 80t


Losen der DG mit dem Euler-Verfahren. Damit von den diskreten Werten im Anschluss die FFT gebildetwerden kann, setzen wir die Anzahl der Zeitschritte als Zweierpotenz> m:=9: N:=2m:

> T:=90.: dt:=T/N:

> t:=0:

> for i from 1 to N

> do

> data1[i]:=[t,U3]: #diskrete Impulsantwort

> data2[i]:=[t,Ue]: #diskretes Eingangssignal

> U1:=U1 + dt*dU1:

> I1 := I1 + dt* dI1:

> U2:= U2 + dt * dU2:

> I2:= I2 + dt * dI2:

> U3:= U3 + dt * dU3:

> t:=t+dt:

> end do:

Graphische Darstellung der Impulsantwort> plot([seq(data1[n],n=1..N)]);

–0.05

0

0.05

0.1

0.15

0.2

20 40 60 80

(2.) Systemanalyse: Um die FFT anwenden zu konnen, benotigen wir die diskreten Funktionswerte fddes Ausgangs- bzw. Eingangssignals. Die Imaginarteile imd sind jeweils auf Null zu setzen.> for i from 1 to N

> do

> fd1[i] := data1[i][2]: #Ausgangssignal

> imd1[i] := 0.:

> fd2[i] := data2[i][2]: #Eingangssignal

> imd2[i] := 0.:

> end do:

6.7 System-Analyse mit Maple 57

Mit der schnellen Fouriertransformation berechnen wir die DFT sowohl des Eingangs- als auch desAusgangssignals berechnet> FFT(m, fd1, imd1);

> FFT(m, fd2, imd2);

Wir bilden anschließend das Verhaltnis von Ausgangs- zu Eingangsspektrum, H(w), und zeichnen denBetrag dieses Verhaltnisses als Funktion von w :> for i from 1 to N

> do H[i]:= evalc( (fd1[i]+I*imd1[i]) / (fd2[i]+I*imd2[i])):

> end do:

> for i from 1 to N

> do daten[i]:= [2*3.14/T*(i-1), abs(H[i])]

> end do:

> plot([seq(daten[i], i=1..N/2)], labels=[‘w‘,“], color=red, thickness=2);

0

0.1

0.2

0.3

0.4

0.5

2 4 6 8 10 12 14 16 18w

Diskussion: Der Graph der Ubertragungsfunktion hat qualitativ den gleichen Verlauf wie der analytischuber komplexe Rechnung berechnete mit der gleichen Grenzfrequenz wg = 1.4.

Aufgabe: Fuhren Sie Simulationen durch, bei denen Sie die Impulsbreite t0 variieren.

Date post:	17-Sep-2018
Category:	Documents
Upload:	truongcong
View:	216 times
Download:	0 times

Numerik von Differenzialgleichungen - home.hs …weth0002/veranstaltungen/details/docs/... ·...

Documents