Wurzelortkurve eines MSV

12.06.2009, 17:55

tigerbine

Wurzelortkurve eines MSV

Gegeben ist folgendes MSV:

$\begin{eqnarray*} \eta_{k+1}:= \eta_k + h [(1-\alpha) f_k + \alpha f_{k+1}] \end{eqnarray*}$

Zu bestimmen ist die Wurzelortkurve in Abhängigkeit von $\begin{eqnarray*} \alpha \end{eqnarray*}$ . Dazu bräuchte ich die charakteristischen Polynome des MSV. Wie bekomme ich die denn? Lese ich nun einfach ab?

$\begin{eqnarray*} a_1=1, a_0 = -1, b_0 = \alpha, b_1 = (1-\alpha) \end{eqnarray*}$

Dann sind

$\begin{eqnarray*} \psi(\lambda) = a_0 \lambda ^0 + a_1 \lambda^1 = -1 + \lambda \end{eqnarray*}$

$\begin{eqnarray*} \chi(\lambda) = b_0 \lambda ^0 + b_1 \lambda^1 = \alpha + (1-\alpha)\lambda \end{eqnarray*}$

Nun steht als Definition der Wurzelortkurve folgendes da:

$\begin{eqnarray*} C=\{\psi(\exp(i\phi)) / \chi(\exp(i\phi))~|~\phi \in [0,2\pi]\} \end{eqnarray*}$

Was bedeutet das denn? Ich setzte den Einheitskreis in die Funktionen ein? Soll ich dann die Funktionswerte gegeneinander auftragen oder dividieren. Sprich was soll "/" bedeuten?

12.06.2009, 21:13

system-agent

Auf diesen Beitrag antworten »

Angenommen du hast das lineare Mehrschrittverfahren
$\begin{eqnarray*} \sum\limits_{i=0}^k\alpha_iy_{n+i}=h\cdot\sum\limits_{i=0}^k\beta_if_{n+i} \end{eqnarray*}$
wobei $\begin{eqnarray*} y_{n} \end{eqnarray*}$ die numerische Lösung am Knotenpunkt $\begin{eqnarray*} n \end{eqnarray*}$ ist und $\begin{eqnarray*} f_{n}:=f(x_n,y_n) \end{eqnarray*}$ .

Dann definiert man die erzeugenden Polynome wie folgt:
$\begin{eqnarray*} \rho(\zeta):=\sum\limits_{i=0}^k\alpha_i\zeta^i \end{eqnarray*}$
$\begin{eqnarray*} \sigma(\zeta):=\sum\limits_{i=0}^k\beta_i\zeta^i \end{eqnarray*}$

[tut mir leid, ich habe immer andere Bezeichnungen...]

Wendet man das MSV auf die Testgleichung $\begin{eqnarray*} y'=\lambda y \end{eqnarray*}$ an, dann bekommt man mit $\begin{eqnarray*} \mu:=h\lambda \end{eqnarray*}$
$\begin{eqnarray*} 0=\sum\limits_{i=0}^k(\alpha_i-\mu\beta_i)y_{n+1} \end{eqnarray*}$
Setze $\begin{eqnarray*} y_n=\zeta^n \end{eqnarray*}$ in dieser Differenzengleichung, bekommt man (für $\begin{eqnarray*} \zeta\neq0 \end{eqnarray*}$ , denn sonst bekommt man triviale Lösungen):
$\begin{eqnarray*} 0=\sum\limits_{i=0}^k(\alpha_i-\mu\beta_i)\zeta^i=\rho(\zeta)-\mu\sigma(\zeta) \end{eqnarray*}$

Als Lösung der Differenzengleichung bekommt man die numerische Lösung und die ist gegeben durch
$\begin{eqnarray*} y_n=c_1(\mu)\zeta_1^n(\mu)+...+c_k(\mu)\zeta_k^n(\mu) \end{eqnarray*}$
wobei die $\begin{eqnarray*} c_i \end{eqnarray*}$ Konstanten sind [für jedes $\begin{eqnarray*} \mu \end{eqnarray*}$ eine] und die $\begin{eqnarray*} \zeta_i \end{eqnarray*}$ die Nullstellen von $\begin{eqnarray*} \rho(\zeta)-\mu\sigma(\zeta) \end{eqnarray*}$ [auch abhängig von $\begin{eqnarray*} \mu \end{eqnarray*}$ ].

Das Stabilitätsgebiet ist dann definiert als:
$\begin{eqnarray*} S:=\{\mu\in\mathbb C\quad:\quad\text{alle Nullstellen von }\rho(\zeta)-\mu\sigma(\zeta)\text{ erfüllen }|\zeta_j(\mu)|\leq1\text{ bzw }|\zeta_j(\mu)|<1\text{ falls mehrfache Nullstelle}\} \end{eqnarray*}$

Das heisst in $\begin{eqnarray*} S \end{eqnarray*}$ sind alle Nullstellen kleiner als 1 [im Betrag]. Betrachte nun den Rand von $\begin{eqnarray*} S \end{eqnarray*}$ . Falls $\begin{eqnarray*} \mu\in\partial S \end{eqnarray*}$ , dann gibt es ein $\begin{eqnarray*} \zeta\in\mathbb C \end{eqnarray*}$ mit $\begin{eqnarray*} |\zeta|=1 \end{eqnarray*}$ so, dass $\begin{eqnarray*} \rho(\zeta)-\mu\sigma(\zeta)=0 \end{eqnarray*}$ .
Das bedeutet, dass es ein $\begin{eqnarray*} \theta\in[0,2\pi] \end{eqnarray*}$ gibt so, dass
$\begin{eqnarray*} e^{i\theta}=\zeta \end{eqnarray*}$ und damit
$\begin{eqnarray*} \rho(e^{i\theta})-\mu\sigma(e^{i\theta})=0 \end{eqnarray*}$ oder umgeformt
$\begin{eqnarray*} \mu=\frac{\rho(e^{i\theta})}{\sigma(e^{i\theta})} \end{eqnarray*}$

Damit [und der Tatsache dass die Nullstellen stetig von $\begin{eqnarray*} \rho(\zeta)-\mu\sigma(\zeta) \end{eqnarray*}$ stetig von $\begin{eqnarray*} \mu \end{eqnarray*}$ abhängen] bekommt man: $\begin{eqnarray*} \partial S\subset\Gamma \end{eqnarray*}$ wobei $\begin{eqnarray*} \Gamma \end{eqnarray*}$ die Spur des Weges
$\begin{eqnarray*} \Gamma:[0,2\pi]\to\mathbb C \end{eqnarray*}$
$\begin{eqnarray*} t\mapsto\frac{\rho(e^{it})}{\sigma(e^{it})} \end{eqnarray*}$

Das ist also die Idee dafür sich eine Hilfe zu verschaffen um das Stabilitätsgebiet zu zeichnen. $\begin{eqnarray*} \Gamma \end{eqnarray*}$ ist ein geschlossener Weg in $\begin{eqnarray*} \mathbb C \end{eqnarray*}$ .
Dieser Weg wird manchmal auch Root Locus Curve genannt.

12.06.2009, 21:49

tigerbine

Auf diesen Beitrag antworten »

Danke für den langen Beitrag. Stimmen den meine Koeffizienten oder ist das falsch was ich gemacht habe?

12.06.2009, 22:41

system-agent

Auf diesen Beitrag antworten »

In dem Fall stimmen deine Koeffizienten $\begin{eqnarray*} a_0 \end{eqnarray*}$ und $\begin{eqnarray*} a_1 \end{eqnarray*}$ aber für die Anderen müsste es $\begin{eqnarray*} b_0=1-\alpha \end{eqnarray*}$ und $\begin{eqnarray*} b_1=\alpha \end{eqnarray*}$ sein.

12.06.2009, 23:35

tigerbine

Auf diesen Beitrag antworten »

Stimmt. (und Mist) Augenzwinkern

Danke, so kann ich mir einige dumme Rechnungen sparen.

Als Bemerkung steht in der Aufgabe, dass Spezialfälle z.B. explizites und Implizites Eulerverfahren wären. Wie würden denn dann die alphas lauten?

expizites Euler:
$\begin{eqnarray*} \eta_{k+1} = \eta_k +h f(t_k,\eta_k) \end{eqnarray*}$

Also $\begin{eqnarray*} \alpha = 0 \end{eqnarray*}$ . Damit folgt $\begin{eqnarray*} a_1=1, a_0 = -1, b_0 = 1, b_1 = 0 \end{eqnarray*}$ .

$\begin{eqnarray*} \psi(\lambda) = a_0 \lambda ^0 + a_1 \lambda^1 = -1 + \lambda \end{eqnarray*}$

$\begin{eqnarray*} \chi(\lambda) = b_0 \lambda ^0 + b_1 \lambda^1 = 1 \end{eqnarray*}$

Dann ist

$\begin{eqnarray*} \frac{\psi(\exp(i\phi))}{\chi(\exp(i\phi))} = \exp(i\phi)-1 \end{eqnarray*}$

Das ist ein Einheitskreis um (-1,0). Genauso wie es sein sollte.

implizites Euler:
$\begin{eqnarray*} \eta_{k+1}=\eta_k+hf(t_{k+1},\eta_{k+1}) \end{eqnarray*}$

Also $\begin{eqnarray*} \alpha = 1 \end{eqnarray*}$ . Damit folgt $\begin{eqnarray*} a_1=1, a_0 = -1, b_0 = 0, b_1 = 1 \end{eqnarray*}$ .

$\begin{eqnarray*} \psi(\lambda) = a_0 \lambda ^0 + a_1 \lambda^1 = -1 + \lambda \end{eqnarray*}$

$\begin{eqnarray*} \chi(\lambda) = b_0 \lambda ^0 + b_1 \lambda^1 = \lambda \end{eqnarray*}$

Dann ist

$\begin{eqnarray*} \frac{\psi(\exp(i\phi))}{\chi(\exp(i\phi))} =1-\frac{1}{\exp(i\phi)} = 1-\exp(-i\phi) \end{eqnarray*}$

Wie erkläre ich nun den Einheitskreis um (1,0). Das "-" in exp ändert ja nur die Drehrichtung. Das - davor spiegelt an der x-Achse. Das macht nun dem Einheitskreis nicht viel. Dann verschiebe ich ihn eben noch um +1. Leider ist das ja nur geometrisch begründet... Also die Bildmenge ist gleich mit der von $\begin{eqnarray*} 1+\exp(i\phi) \end{eqnarray*}$

Hast du einen algebraischen Vorschlag?

Allgemein dann:

$\begin{eqnarray*} \frac{\psi(\exp(i\phi))}{\chi(\exp(i\phi))} =\frac{\exp(i\phi)-1}{1-\alpha + \alpha\exp(i\phi)} \end{eqnarray*}$

Kann man das noch geschickt vereinfachen? Nächste Frage ist, für welche $\begin{eqnarray*} \alpha \end{eqnarray*}$ das ganze A-Stabil ist. Aber das ist doch, zumindest wie hier angegeben komisch. Mit der Wurzelortskurve bekomme ich doch den Rand des Stabilitätsgebiets. Woher weiß man denn, ob es "drinnen" oder "draussen" stabil ist. Die beiden Euler wären da doch ein gutes Beispiel...

Hast du da einen Tipp? Du hast oben was von Nullstellen geschrieben. Das ist in der Anleitung hier gar nicht angegeben. Geht dass denn nur mit dem Rand... verwirrt

13.06.2009, 02:00

system-agent

Auf diesen Beitrag antworten »

Zitat:

Original von tigerbine
expizites Euler:
$\begin{eqnarray*} \eta_{k+1} = \eta_k +h f(t_k,\eta_k) \end{eqnarray*}$

Also $\begin{eqnarray*} \alpha = 0 \end{eqnarray*}$ . Damit folgt $\begin{eqnarray*} a_1=1, a_0 = -1, b_0 = 1, b_1 = 0 \end{eqnarray*}$ .

$\begin{eqnarray*} \psi(\lambda) = a_0 \lambda ^0 + a_1 \lambda^1 = -1 + \lambda \end{eqnarray*}$

$\begin{eqnarray*} \chi(\lambda) = b_0 \lambda ^0 + b_1 \lambda^1 = 1 \end{eqnarray*}$

Dann ist

$\begin{eqnarray*} \frac{\psi(\exp(i\phi))}{\chi(\exp(i\phi))} = \exp(i\phi)-1 \end{eqnarray*}$

Das ist ein Einheitskreis um (-1,0). Genauso wie es sein sollte.

Genau, es ist wie es sein sollte und deshalb auch richtig Freude

Zitat:

Original von tigerbine
implizites Euler:
$\begin{eqnarray*} \eta_{k+1}=\eta_k+hf(t_{k+1},\eta_{k+1}) \end{eqnarray*}$

Also $\begin{eqnarray*} \alpha = 1 \end{eqnarray*}$ . Damit folgt $\begin{eqnarray*} a_1=1, a_0 = -1, b_0 = 0, b_1 = 1 \end{eqnarray*}$ .

$\begin{eqnarray*} \psi(\lambda) = a_0 \lambda ^0 + a_1 \lambda^1 = -1 + \lambda \end{eqnarray*}$

$\begin{eqnarray*} \chi(\lambda) = b_0 \lambda ^0 + b_1 \lambda^1 = \lambda \end{eqnarray*}$

Dann ist

$\begin{eqnarray*} \frac{\psi(\exp(i\phi))}{\chi(\exp(i\phi))} =1-\frac{1}{\exp(i\phi)} = 1-\exp(-i\phi) \end{eqnarray*}$

Wie erkläre ich nun den Einheitskreis um (1,0). Das "-" in exp ändert ja nur die Drehrichtung. Das - davor spiegelt an der x-Achse. Das macht nun dem Einheitskreis nicht viel. Dann verschiebe ich ihn eben noch um +1. Leider ist das ja nur geometrisch begründet... Also die Bildmenge ist gleich mit der von $\begin{eqnarray*} 1+\exp(i\phi) \end{eqnarray*}$

Hast du einen algebraischen Vorschlag?

Ich würde das so belassen. Wems nicht reicht, der soll sich die von dir genannten Abbildung explizit hinschreiben und nachrechnen Augenzwinkern

.

Zunächst gibt die Root Locus Curve wie schon geschrieben nur Auskunft, wo der Rand des Stabilitätsgebietes ist [wie gesagt, er ist Teilmenge der Spur dieser Kurve].
Da die Kurve geschlossen ist [was eine Aufgabe ist zu sehen], umschliesst sie einige Gebiete und man muss eben sehen durch ausprobieren, welches davon zu $\begin{eqnarray*} S \end{eqnarray*}$ gehört. Zum Beispiel bei explizit Euler, man kann den Punkt $\begin{eqnarray*} \mu=1 \end{eqnarray*}$ nehmen. Dann schaut man die Nullstellen von
$\begin{eqnarray*} \rho(\zeta)-\mu\sigma(\zeta)=\rho(\zeta)-\sigma(\zeta) \end{eqnarray*}$
an und schaut die Beträge der Nullstellen an. Sind alle Nullstellen betragsmässig $\begin{eqnarray*} \leq1 \end{eqnarray*}$ und $\begin{eqnarray*} <1 \end{eqnarray*}$ für mehrfache Nullstellen, dann folgt
$\begin{eqnarray*} 1\in S \end{eqnarray*}$ und damit die ganze Kreisscheibe [denn die Kurve gibt schliesslich den Rand...].
Falls $\begin{eqnarray*} 1\notin S \end{eqnarray*}$ gelten würde, dann wüsste man, dass $\begin{eqnarray*} S \end{eqnarray*}$ eben das Komplement der Kreisscheibe wäre [und dann wäre es auch A-Stabil, da dann $\begin{eqnarray*} \mathbb C^- \end{eqnarray*}$ in $\begin{eqnarray*} S \end{eqnarray*}$ wäre; das Argument kannst du für implizit Euler nutzen].

Im Allgemeinen Fall:
Meine erste Idee wäre es zu untersuchen, wann die Root Locus Curve nicht mehr durch $\begin{eqnarray*} \mathbb C^- \end{eqnarray*}$ geht. Das heisst dann zu schauen, wann der Realteil davon positiv ist.

13.06.2009, 10:26

tigerbine

Auf diesen Beitrag antworten »

Ok. Die beiden Euler muss ich nicht bearbeiten. Die waren so als Kontrolle für meine Ideen gedacht. Ich hab nun mal versucht, maple zu bitten, mir einen Plot zu machen... Sehr seltsam... Also ich bekomme eigentlich immer Kreise heraus. Nur bei alpha=0.5 kommt ein seltsames Gebilde zustande. Bis 0.5 liegen die in der neg. Halbebene, dann in der positiven.

Hast du ähnliches raus?

Ich werde dann man schauen, ob ich in dem Buch den Satz mit den "Nullstellen" finde. Dann verwende ich den einfach und verweise darauf.

[attach]10760[/attach][attach]10761[/attach][attach]10762[/attach]

14.06.2009, 22:03

tigerbine

Auf diesen Beitrag antworten »

Jemand einen Tipp?

15.06.2009, 00:30

tigerbine

Auf diesen Beitrag antworten »

Also ich lasse das jetzt bei der Darstellung mit exp. Keine Ahnung, wie man auf die direkte Kreisdarstellung kommt. Meine Rechnungen sind hier zwar richtig (Probepunkte), aber es wird nur häßlicher....

Die Frage nach der Stabilität würde ich, da es sich entweder um explizites Euler oder ein implizites ESV handelt mit der Stabilitätsfunktion beantworten.

$\begin{eqnarray*} R(z)=\frac{1+(1-\alpha)z}{1-\alpha z} \end{eqnarray*}$

Dann mit Argumentation wie hier : A-Stabile Verfahren komme ich dann darauf, dass für $\begin{eqnarray*} \alpha\geq 0.5 \end{eqnarray*}$ das ganze A-Stabil ist. [implizites Euler a=1 spricht zumindest nicht dagegen] Der Weg ist nicht schön, denn ich argumentiere über die Graphen (Parabeln) der Real und Imaginäteile.

15.06.2009, 17:11

tigerbine

Auf diesen Beitrag antworten »

Feedback

Zitat:

Original von tigerbine

$\begin{eqnarray*} \frac{\psi(\exp(i\phi))}{\chi(\exp(i\phi))} =\frac{\exp(i\phi)-1}{1-\alpha + \alpha\exp(i\phi)} \end{eqnarray*}$

Das reicht als explizite Darstellung aus. Auf die Nachfrage, wie man das mit den Kreisen daraus sehen würde, kam die Antwort "Durch plotten"... ROFL

Naja, eigentlich nicht lustig.

Bei $\begin{eqnarray*} \alpha=0.5 \end{eqnarray*}$ ist es die y-Achse. Da hat maple wohl Probleme.

Neue Frage »

Antworten »

Wurzelortkurve eines MSV

Verwandte Themen