Dgl und Numerik Übung 3 mit Lösung PDF

Preview

Summary

Dgl und Numerik für den Maschinenbau von Michael Karow....

Description

L¨ osungen zum 3. DGL-Uebungsblatt Aufgabe 1

4 Punkte

Seien a, ω ∈ R. Bestimmen Sie die allgemeine L¨osung y(t) folgenden DGL. Es handelt sich dabei um inhomogene lineare DGL erster Ordnung mit konstantem Koeffizienten. L¨ osungsmethode ist die Variation der Konstanten. ¨ (Ua) (Ha)

y ′ = a y + t2 y′ = a y + (t + 1)

¨ ( Ub) (Hb)

y′ = a y + sin(ω t) y′ = a y + cos(ω t).

L¨ osungen zu Aufgabe 1 ¨ y(t) = eatc(t) mit c′ (t) = e−at t2 . Somit (Ua) Z c(t) = e−at t2 dt (2mal partiell integrieren um t2 zu eliminieren) Z −at e e−at 2 2t dt = t − −a −a Z e−at 2 2 = − t + e−at t dt a a  Z −at  e e−at 2 2 e−at = − t− t + a −a −a a  2  t 2t 2 = −e−at + 2 + 3 + c. a a a Somit: y(t) =



2t 2 t2 + 2+ 3 a a a



+ eatc.

¨ y(t) = eatc(t) mit c′ (t) = e−at sin(ω t). Zur Berechnung von c(t) wird wieder partiell integriert. (Ub) Diesmal entsteht dabei eine Gleichung f¨ ur c(t), die man nach c(t) umstellt. Ein alternativer L¨ osungsweg (ohne Integration) wird unten angegeben.

1

c(t) = = = = =

Z

e−at sin(ω t) dt (2mal partiell integrieren) Z e−at e−at cos(ω t) ω dt sin(ω t) − −a −a Z e−at ω − e−at cos(ω t) dt sin(ω t) + a a   Z −at e e−at ω e−at cos(ω t) − − sin(ω t) (−ω) dt sin(ω t) + −a −a a a 2 −at −at ω ωe e cos(ωt) − 2 (c(t) + c0 ) sin(ω t) − − a2 a a

⇒ ω2 e−at ω e−at a2 + ω 2 cos(ωt) − 2 c0 c(t) = − sin(ω t) − 2 2 a a a a ⇒   a ω −at c(t) = −e cos(ω t) + c1 sin(ω t) + 2 a + ω2 a2 + ω 2 Somit ist die allgemeine L¨osung der DGL   a ω at y(t) = e c(t) = − cos(ω t) + eatc1 , sin(ω t) + 2 a + ω2 a2 + ω 2

c1 beliebig.

Alternativer L¨ osungsweg: Weil sin(ωt) der Imagin¨arteil von eωt ist und ausserdem a reell ist, kann man folgendermaßen vorgehen: 1) Finde eine partikul¨are L¨osung von yc′ (t) = a yc (t) + eiωt.

(∗)

2) Nimm von der L¨osung yc den Imagin¨arteil: yp (t) = Im(yc (t)). spr¨ unglichen DGL ist dann y(t) = yp (t) + eatc.

Die allgemeine L¨osung der ur-

Um eine partikul¨are L¨ osung von (∗) zu finden, machen wir einen Ansatz vom Typ der rechten Seite: yc (t) = A eiωt. Die Konstante A wird durch Einsetzen in die DGL bestimmt: y′c (t) = A iω eiωt a yc (t) + eiωt = (a A + 1) eiωt

Die rechten Seiten sind genau dann gleich, wenn A iω = a A + 1, also A=

1 −iω − a −iω − a . = = 2 ω + a2 iω − a (iω − a)(−iω − a) 2

Somit:      −iω − a −iω − a −iω − a iω t iω t e Re(eiω t) Im(e ) + Im = Re yp (t) = Im ω 2 + a2 ω 2 + a2 ω 2 + a2 a ω cos(ωt) = − 2 sin(ω t) − 2 2 ω + a2 ω +a Hausaufgaben gehen analog. 

Aufgabe 2 Bestimmen Sie ein reelles Fundamentalsystem zu folgenden DGL. Es handelt sich hierbei um homogene lineare DGL h¨oherer Ordnung mit konstanten Koeffizienten. ¨ ( Ua) ¨ ( Ub) ¨ ( Uc) ¨ ( Ud) (Ha) (Hb) (Hc)

y′′ − 8y′ + γ y = 0, γ ∈ {15, 16, 25} y(4) − 8y(3) + 24y + 8y′ − 25 = 0 y(3) + 6y′′ + 12y + 8y = 0 y(4) = 0 y′′ − 6y′ + γ y = 0, γ ∈ {8, 9, 25} y(4) + 4 y(3) + 6y′′ + 4y′ + y = 0 y(4) + 8y′′ + 16 y = 0

Hinweis: Das charakteristische Polynom von (Hb) ist (λ + 1)4 . L¨ osungen zu Aufgabe 2 Aufgabe ¨ ( Ua) ¨ ( Ua) ¨ ( Ua) ¨ ( Ub) ¨ c) (U ¨ d) (U (Ha) (Ha) (Ha) (Hb) (Hc)

DGL y′′ − 8y′ + 15 = 0 y′′ − 8y′ + 16 = 0 y′′ − 8y′ + 25 = 0 (4) y − 8y(3) + 24y + 8y′ − 25 = 0 y(3) + 6y′′ + 12y + 8 = 0 y(4) = 0 ′′ y − 6y′ + 8 = 0 y′′ − 6y′ + 9 = 0 y′′ − 6y′ + 25 = 0 (4) y + 4 y(3) + 6y′′ + 4y′ + y = 0 y(4) + 8y′′ + 16 = 0

char. Polynom λ2 − 8λ + 15 = (λ − 3)(λ − 5) λ2 − 8λ + 16 = (λ − 4)2 2 λ − 8λ + 25 = (λ − (4 + 3i))(λ − (4 − 3i)) (λ2 − 8λ + 25)(λ2 − 1) (λ + 2)3 λ4 (λ − 2)(λ − 4) (λ − 3)2 (λ − (3 + 4i))(λ − (3 − 4i)) (λ + 2)4 (λ2 + 4)2

Fundamentalsystem e3t , e5t e4t , t e4t e4t sin(3t), e4t cos(3t) e4t sin(3t), e4t cos(3t), et , e−t e−2t , t e−2t , t2 e−2t 1, t, t2 , t3 e2t , e4t e3t , t e3t e3t sin(4t), e3t cos(4t) e−t , e−t t, e−t t2 , e−t t3 sin(2t), cos(2t), t sin(2t), t cos(2t)

Aufgabe 3 Theorie zu Programmieraufgabe. In Anwendungen kommen oft DGL der Form M y′′ + Dy′ + Sy = b(t) vor. Dabei sind M, D, S symmetrische reelle Matrizen gleicher Dimension n (Massenmatrix, D¨ ampfungsmatrix, Steifigkeitsmatrix). M und S sind positiv definit, D ist positiv semidefinit. 1 Im d¨ ampfungsfreien Fall ist D = 0. y = y(t) und b(t) sind vektorwertige Funktionen. In dieser Aufgabe betrachten wir die homogene Gleichung M y′′ + Dy′ + Sy = 0 mit zugeh¨ origem Matrixpolynom 1

(∗)

P (λ) = M λ2 + Dλ + S.

Siehe Foliensammlung zur positiven Definitheit. Skalarer Fall ist das einfache ged¨ a mpfte Federpendel.

3

Ein (m¨ oglicherweise komplexer) Vektor v 6= 0 heisst Eigenvektor von P (·) zum Eigenwert λ ∈ C, wenn P (λ)v = 0. Zeigen Sie (a) F¨ ur jedes Eigenpaar (v, λ) mit λ = α + i ω, α, ω ∈ R, sind die Funktionen y(t) = yR (t) = yI (t) =

eλ t v, Re(y(t)) = eα t (cos(ω t) Re(v) − sin(ω t) Im(v)),

Im(y(t)) = eα t (cos(ω t) Im(v) + sin(ω t) Re(v))

L¨ osungen der DGL. (b) Die Eigenwerte von P sind die Nullstellen des charakteristischen Polynoms χ(λ) = det(P (λ)) = det(M λ2 + Dλ + S). (c) Aufgrund der Definitheit der Matrizen haben alle Eigenwerte nichtpositiven Realteil (Re(λ) ≤ 0). (d) Sei D = 0, also P (λ) = λ2 M + S. Sei A := M −1 S. Dann gilt folgendes.

(d1) Eigenwerte von P sind rein imagin¨ ar: λ = iω, ω ∈ R. Eigenvektoren k¨ onnen reell gew¨ ahlt werden.

(d2) Die Eigenwertgleichung P (λ)v = 0 ist a¨quivalent zu Standard-Eigenwertgleichung f¨ ur die Matrix A: (M −1 S) v = (−λ2 ) v = ω 2 v. |{z} | {z } =:µ

A

(d3) Konsequenz aus (d2): Sei Av = µ v und v reell. Dann ist µ > 0. Setze ω = Funktionen y(t) = eiωt v, yR (t) = cos(ω t)v, yI (t) = sin(ω t)v

√ µ. Dann sind die

L¨ osungen der d¨ ampfungsfreien DGL. Linearkombinationen solcher L¨ osungen zu verschiedenen Eigenwerten/Eigenvektoren sind nat¨ urlich auch wieder L¨ osungen (Superpositionsprinzip). L¨ osungen zu Aufgabe 3

(a) y in DGL einsetzen. Re(·), Im(·) auf DGL anwenden. (b) Die Gleichung P (λ)v = 0 hat genau dann eine von 0 verschiedene L¨osung, wenn det(P (λ)) = 0. Das ist ein Standard-Resultat der linearen Algebra. ⊤ (c) F¨ ur v = vR + i vI , mit reellen vR , vI definiere v∗ := v¯⊤ = v⊤ R − i v I (konjugiert-transponierter Vektor). F¨ ur symmetrische positv definite relle Matrix N und v 6= 0: ⊤ ⊤ ⊤ v∗ N v = (v⊤ R − i vI )N (vR + ivi ) = vR N vR + v I N vI > 0.

Nun: P (λ)v = 0

⇒

0 = v∗ P (λ)v = v∗ M vλ2 + v∗ Dv λ + v∗ Sv.

Das ist eine quadratische Gleichung f¨ ur λ mit positiven Koeffizienten. Mit pq-Formel: s  ∗ 2 v Dv v∗ Dv v∗ Sv λ=− ∗ − ∗ ± ∗ 2 v Mv v Mv 2 v Mv Wenn Term unter der Wurzel negativ ist, dann ist der Term vor der Wurzel der Realteil. Dieser ist ≤ 0. Wenn Term unter Wurzel nicht negativ, dann ist die Wurzel immer noch ≤ als der Betrag der Zahl vor der Wurzel.

4

(d1) Dass Eigenwerte f¨ ur D = 0 imagin¨ ar sind, folgt aus der pq-Formel in der L¨osung von (c). (d2) Mit v = vR + ivI 6= 0, λ = iω , (M (iω )2 + S)(vR + ivI ) = 0

(−M ω 2 + S)vR = 0 und (−M ω 2 + S )vI = 0

⇔

Wegen v 6= 0 ist vR = 6 0 oder vI = 6 0.

d3,d4) sind hoffentlich klar.

Programmieraufgabe. Wir betrachten eine mit Federn verbundene Pendelkette wie im Bild.

y1

y

y

2

y4

3

Das Ziel dieser Aufgabe ist, Filme zu generieren, welche die Pendelbewegung zeigen. Seien mk , sk , yk die Massen, Federsteifigkeiten und Auslenkungswinkel. Sei ℓ die Pendell¨ ange und sei g die Fallbeschleunigung. Bei kleinen Auslenkungswinkeln2 gelten n¨ aherungsweise die DGL m1 ℓy′′1 (t) = m2 ℓy′′1 (t) =

−m1 g y1 (t) + s1 ℓ (y2 (t) − y1 (t)),

m3 ℓy′′3 (t) m4 ℓ y4′′(t)

−m3 g y3 (t) − s2 ℓ (y3 (t) − y2 (t)) + s3 ℓ (y4 (t) − y3 (t)),

= =

−m2 g y2 (t) − s1 ℓ (y2 (t) − y1 (t)) + s2 ℓ (y3 (t) − y2 (t)),

−m4 g y3 (t) − s3 ℓ (y4 (t) − y3 (t)).

Wenn man hier alle Terme auf die linke Seite bringt, dann lauten die DGL in Matrix-Vektor-Schreibweise:     ′′      y1 (t) m1 ℓ 0 y1 (t) m1 g + s 1 ℓ −s 1 ℓ    y2 (t)   0  y2′′ (t)   −s1 ℓ m ℓ m g + ( s + s ) ℓ −s ℓ 2 2 1 2 2     ′′  +        y3 (t)  =  0 .  y3 (t)   m3 ℓ −s 2 ℓ m3 g + (s2 + s3 ) ℓ −s 3 ℓ y4′′ (t) m4 ℓ 0 y4 (t) −s3 ℓ m3 g + s 3 ℓ {z } | {z } | {z } | {z } | M

S

y ′′

Alle nicht ausgef¨ ullten Matrixstellen sind 0. Wie in Aufgabe 3 dargestellt, sind die L¨ osungen     cos(ω1 t) c1 sin(ω1 t) c˜1 4 X cos(ω2 t) c2  sin(ω2 t) c˜2    +  . y(t) = (cos(ω t) vk ck + sin(ω t) vk c˜k ) = [v1 v2 v3 v4]  | {z } cos(ω3 t) c3  sin(ω3 t) c˜3  k=1 V cos(ω4 t) c4 sin(ω4 t) c˜4

y

Dabei sind ck , c˜k beliebig w¨ ahlbare Konstanten. Die vk sind Eigenvektoren von A = M −1 S, und die Frequenzen ωk sind die Wurzeln aus den zugeh¨ origen Eigenwerten: Avk = ω 2k vk . Bei Anfangsgewschingigkeit y′ (0) = 0 enf¨ allt der Sinus-Anteil, und man hat nur   cos(ω1 t) c1  cos(ω2 t) c2  y(t) = [v1 v2 v3 v4]  (∗) | {z }  cos(ω3 t) c3 V cos(ω4 t) c4 2

genauer: wenn yk (t) ≈ sin(yk (t))

5

Setzt man hier nur ein ck 6= 0 , dann bekommt man eine sogenannte Eigenschwingung y(t) = cos(ωk t)vk ck . Alle anderen Schwingungen erh¨ alt man also Linearkombination von Eigenschwingungen (Superpositionsprinzip). Nun zur Aufgabe. Im Aufgabenordner auf ISIS finden Sie das MATLAB-Programm pendelkette.m F¨ ugen Sie dort an angegebener Stelle die Berechnung von V und den ωk ein, und lassen Sie das Programm zun¨ achst f¨ ur die 4 Eigenschwingungen und dann f¨ ur mindestens eine Superposition laufen (mehrere c’s ungleich 0). W¨ ahlen Sie dabei die ck so klein, dass realisitsche Ergebnisse heraus kommen. Benutzen Sie zur Berechnung der Eigenwerte und Eigenvektoren die MATLAB-Funktion eig (Erkl¨ arung: siehe doc eig) Daten: alle Massen=1, alle Steifigkeiten=1, Pendell¨ ange=1, g = 10 Sie d¨ urfen das Programm gerne ver¨ andern oder selbst ein besseres schreiben.

Programmieraufgabe 2. Die Bewegungsgleichung f¨ ur ein einzelnes ebenes mathematisches Pendel mit Auslenkungswinkel y und L¨ ange ℓ ist bekanntlich (Mechanik-Vorlesung, siehe Bild links): y′′ (t) = −

g sin(y) ℓ

y

(∗)

0 y

Die L¨ osungen dieser DGL lassen sich leider nicht mehr mit einer Formel ausdr¨ ucken, in der nur die elementaren Funktionen vorkommen. 3 F¨ ur kleine Auslenkungen y ist aber sin(y) ≈ y (Das erkennt man an der TaylorEntwicklung des Sinus). Die Bewegung ist dann n¨ aherungsweise eine L¨ osung von g y′′ (t) = − y. ℓ

(∗ ∗ ∗)

Diese DGL ist linear und von der gleichen Form wie die DGL f¨ ur ein unged¨ampftes Federpendel (Bild rechts): y′′ (t) = −

s y. m

(m=Masse, s=Steifigkeit)

Die L¨ osungen mit Anfangsgeschwindigkeit y′ (0) = 0 sind p p y(t) = cos(ω t) y(0), ω = g/ℓ bzw. ω = s/m.

(∗ ∗ ∗)

Das St¨ ormer-Verlet-Verfahren ist eine einfache numerische Methode zur L¨osung von DGL der Form y′′ (t) = f (y(t))

(DGL 2. Ordnung)

Das Verfahren basiert auf den Taylor-Entwickungen f¨ ur den Ort y und die Geschwindigkeit v := y′ : y(t + h) =

y(t) + y′ (t) h + |{z}

y′′ (τ1 ) 2 h , 2

v(t)

y (t + h) = y (t) +y′′ (τ2 )h, |{z} | {z } ′

v(t+h)

′

τ1 , τ2 ∈ [t, t + h].

v(t)

Die unbekannten Terme werden wie folgt durch die rechte Seite der DGL ersetzt: y′′ (τ1 ) ≈ y′′ (t) = f (y(t)),

y′′ (τ2 ) ≈

f (y(t)) + f (y(t + h)) . 2

Hier ist ein MATLAB-Programm f¨ ur einen Zeitschritt, um aus bereits berechetem Ort und bereits berechneter Geschwindigkeit (ya,va) (a=alt) auf das neue Paar (yn,vn) zum n¨ achsten Zeitpunkt zu kommen: 3 Mit elementaren Funktionen meine ich Polynome, Exponentialfunktion, trigonometrische Funktionen und ihre Umkehrfunktionen

6

[yn,vn]=stoermerschritt(ya,va,h,f) yn=ya+va*h+f(ya)*h^2/2; vn=va+h*(f(ya)+f(yn))/2; end Nun die Aufgabe. L¨ osen Sie (∗) und (∗∗) mit Hilfe des obigen Verfahrens. Visualisieren Sie die Ergebnisse als Pendel-Film. Visualisieren Sie ebenfalls die exakte L¨ osung (∗ ∗ ∗). Alles in einem Diagramm. Daten: ℓ = 1, g = 10. Schrittweite h und Anfangsauslenkung sind ihre Wahl. Das didaktische Ziel dieser Aufgabe ist u.a. die Unterschiede in den L¨ osungen deutlich zu machen. Zum Plot des Pendels duerfen Sie gern einen Programmteil aus Programmieraufgabe 1 verwenden.

7...