PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : Math - kleine numerische Tricks


richy
28.03.10, 23:15
In diesem Thread mochte ich ein paar einfache numerische Methoden zusammenstellen:


1) Iteratives Erzeugen eines Kreises :
****************************
METHODE A)

Ausgangspunkt sei die Iteration (DZGL) :
z[k+1]:=c*z[k]
Die Loesung lautet
z[k]=z0*c^k= z0*exp( ln(c)*k )
************************

Man sieht an dieser Darstellung nun folgendes :
Ist ln(c) rein imaginaer, so stellt die Loesung einen Kreis in der komplexen Ebene dar :
z[k]=z0*exp( I*s*k )=cos(( I*s*k )+I*sin( I*s*k )

Eine geschlossene Abtastung erhaelt man fuer s=2*Pi/N
Der Hauptwert von ln(c) lautet ln|c|+I(arg(c))
Fuer ln|c|=0 , also |c|=1 wird der Ausdruck rein imaginaer.

Geeinete Faktoren waeren somit
c=c_re+I*c_im=cos(2*Pi/N)+I*sin(2*Pi/N)

Eingesetzt :
z[k+1]:=cos(2*Pi/N)+I*sin(2*Pi/N)*z[k]
******************************

getrennt fuer Real und Imaginaerteil :
z_re[k+1]+I*z_im[k+1]:=(c_re+I*c_im)*(z_re[k]+I*z_im[k])=
c_re*z_re[k]+I*c_re*z_im[k]+I*c_im*z_re[k]-c_im*z_im[k]

Man erhaelt das gekoppelte DZGL Gleichungssystem :
****************************
z_re[k+1]=c_re*z_re[k]-c_im*z_im[k]
z_im[k+1]=c_im*z_re[k]+c_re*z_im[k]
****************************

In Vektorschreibweise
z[k+1]=A*z[k] mit der Matrix A=

*********
c_re , -c_im
c_im , c_re
*********

> restart; with(plots):
> N:=100;
> c_re:=cos(2*Pi/N);
> c_im:=sin(2*Pi/N);
> z_re[1]:=1;
> z_im[1]:=0;
> for k from 1 to N do
> z_re[k+1]:=evalf(c_re*z_re[k]-c_im*z_im[k]);
> z_im[k+1]:=evalf(c_im*z_re[k]+c_re*z_im[k]);
> od:
> druck:=seq([z_re[i],z_im[i]],i=1..N):
> plot ([druck],style=point);

richy
30.03.10, 15:46
Entkoppeln eines partiellen DGL Systems (Charakteristiken Methode)

Notwendige Tools :

DETERMINANTE
***********
http://de.wikipedia.org/wiki/Determinante_%28Mathematik%29
http://upload.wikimedia.org/math/5/5/7/557e8d795884c214008b5b2f486b4902.png

Eigenwert einer Matrix
*****************
Ein Eigenvektor einer Abbildung ist in der linearen Algebra ein vom Nullvektor verschiedener Vektor, dessen Richtung durch die Abbildung nicht verändert wird. Ein Eigenvektor wird also nur gestreckt, und man bezeichnet den Streckungsfaktor als Eigenwert der Abbildung.

http://upload.wikimedia.org/math/7/a/0/7a09c31e1daea4fbae0204f1e11a495f.png
E ist die Einheitsmatrix. Bilden der Determinante fuehrt zu einem Polynom, dass man nach den Eigenwerten lambda aufloest.

Rechtseigenvektor einer Matrix
***********************
http://de.wikipedia.org/wiki/Eigenwertproblem#Berechnung_der_Eigenvektoren
Fuer jeden Eigenwert existiert in der Regel ein Eigenvektor der sich aus folgendem Gleichungssystem bestimmen laesst :
http://upload.wikimedia.org/math/b/c/0/bc011dfb87eee5c295b85c5c4aeb13bf.png

Transponierte einer Matrix
*******************
Beispiel :
http://de.wikipedia.org/wiki/Matrix_%28Mathematik%29#Die_transponierte_Matrix
http://upload.wikimedia.org/math/1/5/7/15767c48d85427c8bf0d5481d2d4d2ea.png
Man vertauscht Zeilen und Spalten

Inverse einer Matrix
***************
http://de.wikipedia.org/wiki/Regul%C3%A4re_Matrix
Beispiel 2x2 Matrix
http://upload.wikimedia.org/math/4/4/9/449b27284b254faa344741a18423e3f6.png

LINKSEIGENVEKTOREN
*****************
.. sind die Rechteigenvektoren der transponierten Matrix

richy
30.03.10, 16:14
ENTKOPPELN EINES allgemeinen 1-D PDE SYSTEMS
**************************************
Das Sytem sei in primitiver Form gegeben :

dw/dt + A(x)*dw/dx + C = 0

Alle Groessen sind Vektoren. A die Uebertragungsmatrix.

TRANSFORMATION AUF DIAGONALGESTALT
********************************
Zunaechst ist es notwendig A(x) auf Diagonalgestalt zu transformieren :
D=S_1*A*S

Dabei ist S eine Matrix, deren Spalten aus den RECHTS-Eigenvektoren der Matrix A bestehen.
S ist die inverse Matrix zu S_1
S_1 besteht aus den transponierten Linkseigenvektoren der Matrix A

Die Elemente der Diagonalmatrix D bestehen dann aus den Eigenwerten der Matrix A

Anwenden auf das PDE System
***********************
Multiplikation von links mit S_1

S_1*dw/dt + S_1*A*dw/x + S_1*C = 0


Trick fuer gemischte Form :
Nun multipliziert man A von rechts Mit der Einheitsmatrix E und drueckt diese aus ueber E=S*S_1

S_1*dw/dt + S_1*A*E*dw/x + S_1*C = 0
S_1*dw/dt + S_1*A*S*S_1*dw/x + S_1*C = 0

Charakteristische Form :
S_1*dw/dt + D*S_1*dw/x + S_1*C = 0
***********************************
Dies ist bereits die entkoppelte PDE in charakteristischen Variablen S_1*w, die sich in einem weiteren Schritt auch in primitiven Variablen umformulieren laesst.
Die Methode arbeitet nur fuer lineare PDE's

richy
30.03.10, 19:09
EDIT : Handrechnung liefert einfachere Ergebnisse als MAPLE

In der Praxis ist oft die Hauptdiagonale mit dem selben Wert belegt. Diesen Fall kann man auch allgemein behandeln :

Beispiel :

Matrix A=
[a11..a12]
[a21..a11]

Eigenwerte :
*********
det (A-L*E) = 0
(a11-L)^2-(a12-a21)=0


L1 := a11-q;
L2 := a11+q
q=(a12*a21)^(1/2)

Linkseigenvektoren :
****************
Sind die Rechtseigenvektoren der transponierten Matrix
Da a12 und a21 in q vertauschbar sind bleiben die Eigenwerte gleich.
Wegen det|(A-L*E)|=0 sind die Gleichungen linear abhaengig. Eine Variable kann als Parameter tau aufgefasst werden. Damit wird eine Gleichung fuer einen beliebigen Wert tau geloest.

Matrix AT (transponiert) =
[a11..a21]
[a12..a11]
****************************
Eigenwert 1 : L1 := a11-q;
(AT-L*E) =
[q..a21]
[a12..-q]

q*x1+a21*tau=0
x1=-a21/q*tau
Ich waehle : tau=x2=-1 =>
Linkseigenvektor 1 = [a21/q, -1]

***************************

****************************
Eigenwert 1 : L2 := a11+q;
(AT-L*E) =
[-q..a21]
[a12..+q]

-q*x1+a21*tau=0
x1=a21/q*tau
Ich waehle : tau=x2=1 =>
Linkseigenvektor 2 = [a21/q, 1]

***************************

Daraus ergibt sich die Linkseigenvektormatrix S_1
[a21/q, -1]
[a21/q, +1]


Die charakteristischen Variablen sind S_1*w :

C1=a21/q*w1 - w2
C2=a21/q*w1 + w2

Die Uebertragungsmatrix ist die Diagonalmatrix D=S_1*A*S
Ihre Diagonale traegt die beiden Eigenwerte L1 und L2

D=
[L1..0]
[0..L2]

SCR
30.03.10, 20:39
Klasse, richy!
(Und ich merke: Ich habe mir da ganz schön was eingebrockt! :D)

richy
30.03.10, 21:11
Hi SCR
Die Rechnung ist noch nicht ganz fertig und loest sich spaeter in Wohlgefallen auf :-)
Man erhaelt eine gemischte Form fuer die Ausgangsvariablen und die Matrix A wird in zwei Matrizen fuer jeden Eigenwert zerlegt. Damit lassen sich sachgemaess harte, weiche und absorbierende Randbedingungen fuer Wellengleichungen formulieren. Ebenso Quellen modellieren. Man koennte auch ein komplettes kleines Wellensimulationsprogramm schreiben.

Fuer den iterativen Kreis wird das Entkoppeln wohl wieder die Ausgangsgleichungen liefern. Vielleicht kann die Methode beim Entscheidungsbaum weiterhelfen. Glaube es zwar fast nicht, aber mal sehen :-)
BTW:
Ich rechne den letzten Thread nochmals von Hand (EDIT)
Insbesonders ist es einfacher zunaechst die Linkseigenvektoren zu berechnen statt die Rechtseigenvektoren.

Gruesse

richy
31.03.10, 00:05
ANWENDUNG ITERATIVER KREIS
************************

Dieser stellt keine PDE dar, aber ein gekoppeltes System.
Viel Neues ist nicht zu erwarten. Insbesonders keine reellen Eigenwerte, denn die Systemfunktion ist eine Drehmatrix.

Iterativer Kreis :
****************************
z_re[k+1]=c_re*z_re[k]-c_im*z_im[k]
z_im[k+1]=c_im*z_re[k]+c_re*z_im[k]
****************************

In Vektorschreibweise
z[k+1]=A*z[k] mit der Matrix A=

*********
c_re , -c_im
c_im , c_re
*********
mit
c_re=cos(2*Pi/N)
c_im=sin(2*Pi/N)


Eigenwerte :
*********
L1 := cos(2*Pi/N)-(-sin(2*Pi/N)^2)^(1/2)
L2 := cos(2*Pi/N)+(-sin(2*Pi/N)^2)^(1/2)

L1:=cos(2*Pi/N)-I*sin(2*Pi/N)
L2:=cos(2*Pi/N)+I*sin(2*Pi/N)


Beide Eigenwerte bilden ueber N jeweils einen Halbkreis in der komplexen Ebene.

Die Transformationsmatrix S_1 gemaess der "Handrechnung" vom letzten Beitrag war :

[a21/q, -1]
[a21/q, +1]

mit
q=Wurzel(a21*a21)
eingesetzt :
q=(-sin(2*Pi/N)^2)^(1/2)

sin(2*Pi/N)^2 ist immer positiv daher die Wurzel stets imaginaer
q=I*sin(2*Pi/N)=c_im=sin(2*Pi/N)=I*a21

q=I*a21

Damit erhalten wir fuer die Transformationsmatrix S_1 :

[I, -1]
[I, +1]


Die charakteristische Variablen lauten dann :

C1=I*w1-w2
C2=I*w1+w2

Fuer C1 ergibt die entkoppelte Gleichung :
C1[k+1]=L1*C1[k]
C1[k+1]=(cos(2*Pi/N)+I*sin(2*Pi/N))*C1[k]

ebenso fuer C2
C2[k+1]=(cos(2*Pi/N)+I*sin(2*Pi/N))*C2[k]

Was bringt uns dieses Ergebnis ?
Fast nichts :-)
Denn wir sind ja schon von einer komplexen, entkoppelten Ausgangsgleichung ausgegangen :
Eingesetzt :
z[k+1]:=cos(2*Pi/N)+I*sin(2*Pi/N)*z[k]
******************************


Das Ergebnis sagt uns daher lediglich :
Bisher haben wir wahrscheinlich richtig gerechnet !
Und den Formalismus koennen wir somit fuer andere Probleme anwenden.

Im naechsten Beitrag stelle ich eine Methode vor, wie man die Matrix A fuer die primitiven Variablen [w1,w2] in zwei Matrizen M1 und M2 zerlegt, die jeweils die Eigenwerte und Eigenvektoren der Matrix beinhalten. Also eine gemischt primitive, charakteristische Form.
So long.

richy
31.03.10, 01:00
Ok, lets go
Wo waren wir in der Theorie stehen geblieben ?
Hier :
Charakteristische Form :
S_1*dw/dt + D*S_1*dw/x + S_1*C = 0
Die Linkseigenvektormatrix S_1 haben wir bereits berechnet.
[a21/q, -1]
[a21/q, +1]
Die Diagonalmatrix D traegt die Eigenwerte.

Was machen wir nun ?
Wir fassen D*S_1*dw/x in einem Vektor M zusammen :
=>
S_1*dw/dt + M + S_1*C = 0 ...

Der Vektor M ergibt sich fuer die spezielle Matrix A zu:
M1=(a11-q)*a21/q*dw1/dx-(a11-q)*dw2/dx
M2=(a11+q)*a21/q*dw1/dx-(a11+q)*dw2/dx

... und transformieren das System wieder auf seine primitiven Variablen zurueck, indem wir es linksseitig mit der zu S_1 inversen Matrix S multiplizieren :
S*S_1*dw/dt + S*M + S*S_1*C = 0 ...

dw/dt + S*M + C = 0

Was ist das fuer ein seltsamer Trick ?
In S*M ist noch die charakteristische Form enthalten, so dass wir die Matrix S*M in ihre charakteristischen Anteile zerlegen koennen.

S laesst sich ueber die Inverse von S_1 bestimmen.
Oder ueber einen Hilfsvektor
d=S*M.
S_1*d=M

Mit der Cramerschen Regel erhaelt man fuer den Hilfsvektor d fuer die spezielle Matrix A :
d1= q/a21/2*(M1+M2)
d1= 1/2*(-M1+M2)

und schliesslich die gemischte Form :

dw1/dt+q/a21/2*( M1 + M2 ) +C1 = 0
dw2/dt+1/2*(-M1 + M2 ) +C2 = 0

mit

M1=(a11-q)*a21/q*dw1/dx-(a11-q)*dw2/dx
M2=(a11+q)*a21/q*dw1/dx+(a11+q)*dw2/dx

wobei d/dt und d/dx durch beliebige Operatoren ersetzt werden koennen.

Anwendungsversuch auf den iterativen Kreis im naechsten Beitrag.

richy
31.03.10, 20:46
Hi
dw/dt + A(x)*dw/dx + C = 0

dw1/dt+q/a21/2*( M1 + M2 ) +C1 = 0
dw2/dt+1/2*(-M1 + M2 ) +C2 = 0

M1=(a11-q)*a21/q*dw1/dx-(a11-q)*dw2/dx
M2=(a11+q)*a21/q*dw1/dx+(a11+q)*dw2/dx


In einer Wellengleichung waeren M1 und M2 jeweils die einlaufenden und Auslaufenden Wellenanteile. (Eine Wellengleichung hat stets zwei Loesungen)
Ueber M1 und M2 lassen sich diese Anteile vorgeben. Bei einer Wellengleichung kann man z.B. eine Ausbreitungsrichtung unterdruecken (Absorbtion am Rand).
Laesst sich das System entkoppeln, so stellt jede Gleichung einen Transportvorgang dar. Ueber eine Koordinatentransformation ueber die man sich in der Transportgeschwindigkeit mitbewegt wird d/dt=0.
Mittels der Charakteristikenmethode lassen sich somit einfache Wellengleichungen auch analytisch loesen. In krummlinigen Koordinaten laesst sich in der Regel das System nicht mehr entkoppeln. Es treten geometrische Quellterme in der Konstanten C auf. Das System bleibt darueber verkoppelt.
Gibt man im iterativen Kreis spezielle Bedingungen fuer M vor, so werden die Variablen selbst komplexwertig und erzeugen einen Kreis. Das ist weniger interessant.

richy
31.03.10, 21:23
Zur Abwechslung etwas einfacheres :

IN DER VERTEILUNG PROGRAMMIERBAHRER ZUFALLSGENERATOR
***********************************************

Die Random Algorithmen von Programmiersprachen erzeugen in der Regel gleichverteilte Zufallszahlen. Wie geht man vor, wenn man einen Zufallszahlengenerator benoetigt, der beliebig verteilte Zufallszahlen erzeugt ? Diese Aufgabestellung ist etwas schwieriger als sie zunaechst erscheint. Man kann sie mit Hilfe des aequivalenten Ereignisses loesen.

Die Ziel-Verteilung wird dabei ueber Abtastwerte im Array a[N] vorgegeben.
Im Beispielprogramm ist die Variable j (nach dem break in der for Schleife) die Variable,
die mit der in a[N] gespeicherten Haufigkeit auftritt.
Ein recht nuetzliches Tool.

Hier der Quellcode in Maple:
**********************
# programmierbarer randomgenerator
restart;
with(linalg):
N:=6;
N_go:=5000; # Anzahl Simulationsschritte
z:=rand(1..1000)/1000: # z=Zufallszahl 0..1
for i from 1 to N+1 do b[i]:=0; od: # testfeld

# Vorgabe Wahscheinlichkeit;
a[1]:=1;a[2]:=2;a[3]:=3;a[4]:=4;a[5]:=5;a[6]:=6;

#Normieren
p:=0:
for i from 1 to N do p:=p+a[i]; od:
for i from 1 to N do a[i]:=a[i]/p; od:
# Integrieren #
a[1]:=a[1]; for i from 2 to N do a[i]:=a[i]+a[i-1]; od;

# SIMULATION. Hier werden die Zufallszahlen erzeugt
for i from 1 to N_go do
s:=z(): # s wird die Zufallsfunktion zugeordnet
for j from 1 to N do
if s<=a[j] then break; fi; # s=Zufallszahl
od;
b[j]:=b[j]+1; # Zaehlen ob der Generator funktioniert
od:

# GRAPHISCHE DARSTELLUNG, yepp er funktioniert :-)
druck2:=seq([(i),(b[i])],i=1..N+1):
plot([[druck2]])

Fuer einige spezielle bekannte Verteilungen existieren besondere spezielle Algorithmen.