Zurück zum Inhaltsverzeichnis

Kapitel 26

Integraltransformationen und Z-Transformation

Maple stellt im Paket inttrans folgende Transformationen zur Verfügung:

> with(inttrans);

[addtable, fourier, fouriercos, fouriersin, hankel,...
[addtable, fourier, fouriercos, fouriersin, hankel,...

Das sind also drei Fouriertransformationen, sowie die Laplace-, Mellin-, Hankel- und Hilberttransformation. Mit 'addtable' und 'savetable' können zu den jeweiligen Transformationen Namensgebungen und Rechenregeln hinzugefügt und abgespeichert werden. Mit ?examples,index finden Sie Beispiele unter 'Integral Transforms'.

Die Fouriertransformation wurde wegen ihrer herausragenden Bedeutung (und als Bindeglied Reihenentwicklung - Integraltransformation) schon im vorangehenden Kapitel behandelt.

In diesem Kapitel steht bei den Integraltransformationen die Laplacetransformation im Vordergrund. Im Zeitalter der Digitalisierung darf aber ihre diskretisierte Verwandte, die Z-Transformation, nicht fehlen. (Die Z-Transformation befindet sich nicht im Paket inttrans.)

>

Laplacetransformation

Bei der Laplacetransformation (kurz L-Transformation) wird vorausgesetzt, dass die Originalfunktion f(t) für t<0 Null ist und für t>=0 stückweise glatt. Die Laplacetransformation ist folgendermaßen definiert:

> restart:with(inttrans):

> F(s):=convert(laplace(f(t), t, s),int);

F(s) := int(f(t)*exp(-t*s),t = 0 .. infinity)

Im Gegensatz zur Fouriertransformation ist der Exponent des Kerns nun nicht mehr rein imaginär, sondern s ist eine komplexe Zahl (mit einem Realteil), es kann also Probleme mit der Konvergenz des Integrals geben. Wenn die Originalfunktion wie exp(a*t) wächst, so ist die Bildfunktion nur für Re(s)>a definiert. Maple kann dem Benutzer die nun erforderlichen Fallunterscheidungen nicht abnehmen. Ein einfaches Beispiel mag dies verdeutlichen:

> f:=t->1;

f := 1

> value(F(s));

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> s

Will now try indefinite integration and then take limits.

limit(-(exp(-t*s)-1)/s,t = infinity)

> assume(s>0):%;

1/s

Bitte behalten Sie dies im Hinterkopf, wenn Maple entweder kein Ergebnis liefert, oder nicht das erwartete. Man hört immer wieder (oder liest es in den Mailinglisten): 'Warum vereinfacht Maple den Term nicht so wie ich mir das vorstelle?'. In den meisten Fällen kann das aber darauf zurückgeführt werden, daß der Benutzer in reellen Zahlen denkt oder sich um die Konvergenz eines Integrals nicht kümmert, während Maple für solche Fälle Stoppstellen eingebaut hat, über die man nur kommt, wenn man die passenden Annahmen macht (der Dialog zwischen CAS und Benutzer läßt sich mit ?infolevel fördern).

Nun ein paar Beispiele zum Gebrauch von 'laplace' und 'invlaplace':

> s:='s':

> laplace(1,t,s);

1/s

> invlaplace(%,s,t);

1

> laplace(t^n,t,s);

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> n+1

Will now try indefinite integration and then take limits.

laplace(t^n,t,s)

So ähnlich hatten wir das doch schon einmal?

> assume(n,posint):

> p:=laplace(t^n,t,s);

p := s^(-n-1)*GAMMA(n+1)

> invlaplace(p,s,t);

t^n

Testen Sie auch andere Annahmen wie z.B. assume(n,negint): und assume(n,nonposint):

Wegen der Voraussetzung, dass die Originalfunktion f(t) für t<0 verschwindet, spielt die Heavisidefunktion im Zusammenhang mit der L-Transformation eine wichtige Rolle

>

> alias(sigma=Heaviside):

> laplace(sigma(t),t,s);

1/s

Das ist die Transformierte von 1, weil Maple davon ausgeht, dass der Benutzer mit 1 meint: 'f(t)=1 für t>=0 und Null für t<0'. Die verschobene Sprungfunktion erfordert aber bereits eine Fallunterscheidung:

> a:='a':laplace(sigma(t-a),t,s);

laplace(sigma(t-a),t,s)

> assume(a>0):%;

exp(-s*a)/s

> invlaplace(%,s,t);

sigma(t-a)

> assume(a<0):%%%;

1/s

> invlaplace(%,s,t);

1

Die Verschiebung nach links geht also bei der Rücktransformation verloren, weil die Originalfunktion laut Annahme immer links von Null verschwindet.

>

Weitere Beispiele zur Verschiebung von Funktionen

So wird eine Funktion korrekt verschoben (Verschiebungssatz):

> f:='f':versch:=laplace(f(t-a)*sigma(t-a),t,s);

versch := laplace(f(t-a)*sigma(t-a),t,s)

> assume(a>0):%;

exp(-s*a)*laplace(f(t),t,s)

> assume(a<0):%%;a:='a':

laplace(f(t-a),t,s)

>

So wird eine Funktion abgeschnitten:

> abschn:=laplace(f(t)*sigma(t-a),t,s);

abschn := laplace(f(t)*sigma(t-a),t,s)

> assume(a>0):%;

exp(-s*a)*laplace(f(t+a),t,s)

> assume(a<0):%%;a:='a':

laplace(f(t),t,s)

>

So wird meistens nicht das erreicht, was man will:

> linears:=laplace(f(t-a),t,s);

linears := laplace(f(t-a),t,s)

> assume(a>0):%;

laplace(f(t-a),t,s)

> assume(a<0):%%;a:='a':

laplace(f(t-a),t,s)

>

Konkretes Beispiel

> f:=x->x^2; # hier können eigene Funktionen getestet werden

f := proc (x) options operator, arrow; x^2 end proc...

> versch;

diff(laplace(sigma(t-a),t,s),`$`(s,2))+2*a*diff(lap...

D.h., es wird zunächst nach dem Linearitätssatz ausgewertet

> assume(a>0): simplify(versch);

2*exp(-s*a)/(s^3)

In Übereinstimmung mit dem Verschiebungssatz (Verschiebung nach rechts)

> assume(a<0):simplify(versch);

(2-2*a*s+a^2*s^2)/(s^3)

In Übereinstimmung mit dem Verschiebungssatz (Verschiebung nach links, oben ist a<0)

> simplify(exp(a*s)*(laplace(f(t),t,s)-int(exp(-s*t)*f(t),t=0..a)));

(s^2*a^2+2*s*a+2)/(s^3)

>

Abschneiden der unverschobenen Funktion rechts vom Ursprung ist keine Verschiebung nach links:

> assume(a>0): simplify(abschn);

exp(-s*a)*(a^2*s^2+2*s*a+2)/(s^3)

und hat links vom Ursprung wegen der Annahme f(t)=0 für t<0 keine Wirkung

> assume(a<0):simplify(abschn);

2*1/(s^3)

>

Eine Verschiebung der Originalfunktion ohne Berücksichtigung von f(t)=0 für t<0 ist nicht im Sinne des Verschiebungssatzes

> assume(a>0): simplify(linears);

(2-2*a*s+a^2*s^2)/(s^3)

was aber bei einer Verschiebung nach links keinen Fehler verursacht

> assume(a<0):simplify(linears);

(2-2*a*s+a^2*s^2)/(s^3)

>

Die zwei wichtigsten Eigenschaften der L-Transformation

1. Die Transformierte von Ableitungen (und Integralen)

> restart:with(inttrans):

> alias(F(s)=laplace(f(t),t,s)):

> laplace( diff(f(t), t), t, s);

s*F(s)-f(0)

> laplace( diff(f(t), t$2), t, s);

s*(s*F(s)-f(0))-D(f)(0)

> laplace(int(f(T), T=0..t), t, s);

F(s)/s

D.h., die L-Transformierte der Ableitung einer Funktion ist s mal die L-Transformierte der Funktion (minus Anfangswert). Damit können Differentialgleichungen algebraisiert werden:

>

> dgl:=diff(f(t), t$2)+ diff(f(t), t)+f(t)=0;

dgl := diff(f(t),`$`(t,2))+diff(f(t),t)+f(t) = 0

Transformation

> alggl:=laplace( dgl ,t,s);

alggl := s*(s*F(s)-f(0))-D(f)(0)+s*F(s)-f(0)+F(s) =...

Lösung im Bildbereich

> bildsol:=solve(alggl,F(s));

bildsol := (s*f(0)+D(f)(0)+f(0))/(s^2+s+1)

Rücktransformation

> sol:=invlaplace(bildsol,s,t);

sol := exp(-1/2*t)*f(0)*cos(1/2*sqrt(3)*t)+1/3*exp(...

Das läuft also ab, wenn man dsolve mit method=laplace aufruft:

> dsolve(dgl,f(t),method=laplace);

f(t) = exp(-1/2*t)*f(0)*cos(1/2*sqrt(3)*t)+1/3*exp(...

Zum Vergleich

> dsolve(dgl,f(t));

f(t) = _C1*exp(-1/2*t)*cos(1/2*sqrt(3)*t)+_C2*exp(-...

Vordergründig besteht der Vorteil der Laplacemethode gegenüber 'dsolve-pur' also darin, dass man die Anfangsbedingungen 'frei Haus' bekommt und sich nicht mit den Integrationskonstanten herumschlagen muß. Der tiefere Sinn der Laplacemethode liegt aber in der Anpassung der Methode an das Problem. Diese Anpassung kann noch weiter verfeinert werden, wenn der Benutzer die Methode nicht über dsolve aufruft, sondern obige Schritte selbst durchführt, um im Bedarfsfall Maple helfen zu können.

>

2. Faltungssatz

> restart:with(inttrans):

Zur besseren Übersicht werden folgende Abkürzungen verwendet:

> alias(L=laplace):alias(F1(s)=laplace(f1(t),t,s)):alias(F2(s)=laplace(f2(t),t,s)):

>

Zur L-Transformation gehört die einseitige Faltung (von 0 bis t)

> faltung:=int(f1(t-tau)*f2(tau),tau=0..t);

faltung := int(f1(t-tau)*f2(tau),tau = 0 .. t)

> bild:=L(faltung,t,s);

bild := F2(s)*F1(s)

In Worten: Die Transformierte einer Faltung ist das Produkt der Transformierten - Algebraisierung!

>

Auch die Umkehrung ist Maple bekannt:

> invlaplace(bild,s,t);

int(f2(_U1)*f1(t-_U1),_U1 = 0 .. t)

>

Anwendung siehe Lösung von DGLn mit der L-Transformation

>

Für Faltungen wie sie bei Messungen vorkommen ist die L-Transformation weniger geeignet, weil die L-Faltung einseitig ist. Beispiele zur zweiseitigen Faltung finden Sie bei der Fouriertransformation.

>

Addtable, savetable

Zu addtable finden Sie einige Beispiele in der Maple-Hilfe: ?examples/addtable

Hier ist noch ein Beispiel wie man neben Namensgebungen und Rechenregeln konkrete Einträge hinzufügt:

> restart:with(inttrans):

> laplace(exp(-x^2),x,s);

1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s)

> invlaplace(%,s,x);

1/2*sqrt(Pi)*invlaplace(exp(1/4*s^2)*erfc(1/2*s),s,...

Die Rücktransformierte ist zunächst nicht bekannt, kann aber mit addtable hinzugefügt werden

> addtable(invlaplace,exp(1/4*s^2)*erfc(1/2*s),2/sqrt(Pi)*exp(-x^2),s,x);

Mit infolevel kann man Maple beim Suchen zuschauen

> infolevel[all]:=5:

> invlaplace(1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s),s,x);

inttrans/invlaplace/lookup:   Entered table lookup (invlaplace)   exp(1/4*s^2)*erfc(1/2*s)

inttrans/invlaplace/lookup:   Looking up   exp(1/4*s^2)*erfc(1/2*s)   in Whittaker table

inttrans/invlaplace/lookup:   Looking up   exp(1/4*s^2)*erfc(1/2*s)   in other table

inttrans/invlaplace/lookup:   Returning from lookup (invlaplace)   2/Pi^(1/2)*exp(-x^2)

exp(-x^2)

Beim zweiten Aufruf wird die Information nicht mehr angezeigt

> invlaplace(1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s),s,x);

exp(-x^2)

weil der Eintrag nun auch in der remember table von invlaplace steht (dort schaut Maple zuerst nach)

> op(4,eval(invlaplace));

TABLE([(1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s), s, x...

>

Bitte beachten Sie, daß der Eintrag ohne den konstanten Vorfaktor gemacht werden muß:

> restart:with(inttrans):

> laplace(exp(-x^2),x,s);

1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s)

> invlaplace(%,s,x);

1/2*sqrt(Pi)*invlaplace(exp(1/4*s^2)*erfc(1/2*s),s,...

> addtable(invlaplace,1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s),exp(-x^2),s,x);

>

> infolevel[all]:=2:

> invlaplace(1/2*sqrt(Pi)*exp(1/4*s^2)*erfc(1/2*s),s,x);infolevel[all]:='infolevel[all]':

inttrans/invlaplace/lookup:   Entered table lookup (invlaplace)   exp(1/4*s^2)*erfc(1/2*s)

simplify:   applying  commonpow  function to expression

simplify:   applying  power  function to expression

combine:   combining with respect to   combine/exp

combine:   combining with respect to   combine/exp

simplify:   applying   simplify/power   function to expression

inttrans/invlaplace/lookup:   Entered table lookup (invlaplace)   exp(1/4*tt^2)

combine:   combining with respect to   combine/exp

simplify:   applying   simplify/power   function to expression

inttrans/invlaplace/lookup:   Entered table lookup (invlaplace)   erfc(1/2*tt)

simplify:   applying   simplify/power   function to expression

inttrans/invlaplace/lookup:   Entered table lookup (invlaplace)   exp(1/4*tt^2)*erfc(1/2*tt)

1/2*sqrt(Pi)*invlaplace(exp(1/4*s^2)*erfc(1/2*s),s,...

>

Die Tabelleneinträge können mit savetable gespeichert und wieder eingelesen werden werden.

>

>

Lösung von Differentialgleichungen mit der Laplacetransformation

Sprungantwort eines RC-Glieds

Das prinzipielle Vorgehen beim Lösen von DGLn wurde schon beschrieben. Hier sind weitere konkrete Beispiele.

Die Sprungantwort eines RC- oder RL-Glieds kann durch eine DGL erster Ordnung beschrieben werden

> restart:with(inttrans): alias(sigma=Heaviside):

> de:=diff(y(t),t)+b*y(t)=u;

de := diff(y(t),t)+b*y(t) = u

> u:=sigma(t-1)-sigma(t-5);

u := sigma(t-1)-sigma(t-5)

Transformsation der DGL

> alias( Y(s)=laplace(y(t),t,s) ):

> laplace(de, t, s);

s*Y(s)-y(0)+b*Y(s) = exp(-s)/s-exp(-5*s)/s

Auflösen nach der Transformierten

> sol:=solve(subs(y(0)=0,%), Y(s));

sol := (exp(-s)-exp(-5*s))/(s*(s+b))

Rücktransformation

> antwort:=invlaplace(sol ,s,t);

antwort := sigma(t-1)/b-sigma(t-1)*exp(-b*(t-1))/b-...

>

> b:=2:

> plot({b*antwort,u}, t=0..9);

[Maple Plot]

>

RLC-Kreis

Ebenfalls eine bekannte DGL, diesmal explizit mit der L-Transformation gelöst

> restart:with(inttrans):

> sys:= u(t)=R*j(t) + L*diff(j(t),t) + 1/C*int(j(t), t):

> de:= diff(sys, t);

de := diff(u(t),t) = R*diff(j(t),t)+L*diff(j(t),`$`...

> alias( U(s)=laplace(u(t), t, s), J(s)=laplace(j(t), t, s)):

Transformation

> bild:=laplace( de, t, s);

bild := s*U(s)-u(0) = R*(s*J(s)-j(0))+L*(s*(s*J(s)-...

Auflösen nach dem transformierten Strom

> Js:=solve(bild, J(s));

Js := C*(s*U(s)-u(0)+R*j(0)+L*s*j(0)+L*D(j)(0))/(R*...

Einsetzen von Anfangsbedingungen

> Js:=subs( u(0)=0, j(0)=0, D(j)(0)=0, Js);

Js := C*s*U(s)/(R*C*s+L*C*s^2+1)

Bereitstellen einer transformierten Spannung (hier im Originalbereich konstant)

> U(s):=laplace(ut, t, s);

U(s) := ut/s

> jt:=invlaplace(Js, s, t);

jt := 2*ut*exp(-1/2*R*t/L)*sin(1/2*sqrt((4*L-R^2*C)...

Werte einsetzen

> R:=12: L:=2: C:=0.02: ut:=200:

> plot({ut, jt*R}, t=0..2);

[Maple Plot]

Sinusförmige Erregerspannung

> R:='R':C:='C':L:='L':U(s):='U(s)':

> ut:=100*sin(2*Pi*50*t):

> us:=laplace(ut, t, s);

us := 10000*Pi/(s^2+10000*Pi^2)

> invlaplace( subs( U(s)=us, Js), s, t):

> jt:=subs(R=10, L=2, C=0.02, %): simplify(%): evalf(%,5);

-.15915*cos(314.16*t)+.25337e-2*sin(314.16*t)+.1591...
-.15915*cos(314.16*t)+.25337e-2*sin(314.16*t)+.1591...

> plot( {jt*10}, t=0..1, numpoints=2000);

[Maple Plot]

Trapezförmige Erregerspannung

> alias(sigma=Heaviside):

> ut:=t-(t-1)*sigma(t-1)-(t-4)*sigma(t-4)+(t-5)*sigma(t-5);

ut := t-(t-1)*sigma(t-1)-(t-4)*sigma(t-4)+(t-5)*sig...

> us:=laplace(ut, t, s);

us := 1/(s^2)-exp(-s)/(s^2)-exp(-4*s)/(s^2)+exp(-5*...

Werte einsetzen

> subs( R=8, L=2, C=0.02,U(s)=us, Js);

.2e-1*s*(1/(s^2)-exp(-s)/(s^2)-exp(-4*s)/(s^2)+exp(...

> invlaplace(%,s,t): jt:=evalf(%,5);

jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...
jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...
jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...
jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...
jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...
jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...
jt := .20000e-1-.87287e-2*exp(-2.*t)*sin(4.5826*t)-...

> plot( {ut, jt*16}, t=0..8);

[Maple Plot]

>

Wärmeleitung (partielle DGL)

Eindimensionale Wärmeleitungsgleichung (homogenes Medium, b ist der Leitfähigkeit proportional):

> restart:with(inttrans):

> pde:=b^2*diff(u(t,x),x$2)-diff(u(t,x),t)=0;

pde := b^2*diff(u(t,x),`$`(x,2))-diff(u(t,x),t) = 0...

Lösungsversuch mit pdsolve

> sol:=pdsolve(pde,u);

sol := `&where`(u(t,x) = _F1(t)*_F2(x),[{diff(_F1(t...

> solb:=pdsolve(pde,u,build);

solb := u(t,x) = _C1*exp(_c[1]*t)*_C2*exp(sqrt(_c[1...

Kann physikalisch nicht sein: _c[1] müsste negativ sein, das würde aber eine Oszillation mit x bedeuten.

L-Transformation

> alias(U=laplace(u(t,x),t,s)):

> bildgl:=laplace(pde,t,s);

bildgl := b^2*diff(U,`$`(x,2))-s*U+u(0,x) = 0

Für die Lösung der DGL benötigt dsolve eine Funktion von x, die Anfangsbedingung wird 0 gesetzt (kann nach der Rücktransformation wieder hergestellt werden).

> de:=subs(U=T(x),u(0,x)=0,bildgl);

de := b^2*diff(T(x),`$`(x,2))-s*T(x) = 0

Transformierte Randbedingungen A(s)

> bildsol:=dsolve({de,T(0)=A0(s),T(l)=A1(s)},T(x));

bildsol := T(x) = (exp(-sqrt(s)*l/b)*A0(s)-A1(s))*e...

Ranbedingung mit l -> infinity

> bildsoll:=limit(bildsol,l=infinity);

bildsoll := T(x) = limit((exp(-sqrt(s)*l/b)*A0(s)-A...

Kann nur ausgewertet werden, wenn über s und b nähere Angaben gemacht werden

> assume(s>0,b>0):bildsolla:=bildsoll;

bildsolla := T(x) = exp(-sqrt(s)*x/b)*A0(s)

Maple benötigt für den Faltungssatz noch die Angabe, dass A0(s) eine transformierte Randbedingung ist

> A0(s):=laplace(a0(t),t,s):

> sol:=invlaplace(rhs(bildsoll),s,t);

sol := int(1/2*a0(_U1)*x*exp(-1/4*x^2/(b^2*(t-_U1))...

Angabe einer Randbedingung im Originalbereich

> subs(a0(_U1)=a0,sol);

int(1/2*a0*x*exp(-1/4*x^2/(b^2*(t-_U1)))/(b*sqrt(Pi...

> slov:=value(%);

slov := limit(a0*erf(1/2*x/(b*sqrt(t-_U1)))-erf(1/2...

Auch hier wieder eine nähere Angabe

> assume(x>0):solv:=%;

solv := a0-erf(1/2*x/(sqrt(t)*b))*a0

Nun müssen die Annahmen wieder entfernt werden, damit mit der Lösung weitergearbeitet werden kann, diesmal mit einem Trick:

> solv:=parse(convert(solv,string));

solv := a0-erf(1/2*x/(sqrt(t)*b))*a0

Konkrete Werte

> b:=0.5:a0:=1:

> plot3d(solv,t=0.001..10,x=-5..5,axes=framed,grid=[10,50],style=patchcontour,contours=20,shading=Z);

[Maple Plot]

>

Weitere Integraltransformationen

Hier werden zur Vollständigkeit alle Transformationen aufgelistet, die in Maple bekannt sind (außer der Laplacetransformation), zusätzlich wird jeweils die Rücktransformation und die Transformierte von diff(f(t),t) angegeben (sofern verfügbar). Beispiele finden Sie in ?examples,index -> Integral Transforms.

Die Fouriertransformation

> restart: with(inttrans):

> convert(fourier(f(t),t,s),int);

int(f(t)*exp(-I*t*s),t = -infinity .. infinity)

und ihre Rücktransformation

> convert(invfourier(f(t),t,s),int);

1/2*int(f(t)*exp(I*t*s),t = -infinity .. infinity)/...

> fourier(diff(f(t),t),t,s);

I*s*fourier(f(t),t,s)

>

Fouriersinus

> convert(fouriersin(f(t),t,s),int);

sqrt(2)*int(f(t)*sin(t*s),t = 0 .. infinity)/(sqrt(...

ist selbstinvertierend

> fouriersin(fouriersin(f(t),t,s),s,t);

f(t)

> fouriersin(diff(f(t),t),t,s);

-s*fouriercos(f(t),t,s)

>

Fouriercosinus

> convert(fouriercos(f(t),t,s),int);

sqrt(2)*int(f(t)*cos(t*s),t = 0 .. infinity)/(sqrt(...

> fouriercos(fouriercos(f(t),t,s),s,t);

f(t)

> fouriercos(diff(f(t),t),t,s);

(-sqrt(2)*f(0)+s*fouriersin(f(t),t,s)*sqrt(Pi))/(sq...

(enthält eine Anfangsbedingung)

>

Mellintransformation

> convert(mellin(f(t),t,s),int);

int(f(t)*t^(s-1),t = 0 .. infinity)

Zur Rücktransformation siehe ?examples,mellin

insbesondere:

"The lesson to be learned from the above is that you must have the correct assumptions on parameters for all integral transforms, and the correct range for the inverse Mellin transform."

>

Hankeltransformation

> convert(hankel(f(t),t,s,v),int);

int(f(t)*sqrt(t*s)*BesselJ(v,t*s),t = 0 .. infinity...

> v:='v':

> hankel(hankel(f(t),t,s,v),s,z,v);

hankel(hankel(f(t),t,s,v),s,z,v)

> assume(v>-1/2,v<>0):%;

f(z)

>

> hankel(diff(f(t),t),t,s,v);

-1/2*s*(hankel(f(t),t,s,v-1)*v+hankel(f(t),t,s,v-1)...

>

Hilberttransformation

> convert(hilbert(f(t),t,s),int);

int(f(t)/(t-s),t = -infinity .. infinity)/Pi

> invhilbert(f(s),s,t);

-hilbert(f(s),s,t)

> hilbert(hilbert(f(t),t,s),s,t);

-f(t)

> hilbert(diff(f(t),t),t,s);

diff(hilbert(f(t),t,s),s)

>

Z-Transformation

Die Z-Transformation steht im Maple-Browser zwar unter Integraltransformationen, ist aber eine diskrete Transformation und kann auch angesprochen werden, ohne das Paket inttrans zu laden.

Einführung

Die Z-Transformation

> restart:

> F(z) = sum(f(n)/z^n,n=0..infinity);

F(z) = sum(f(n)/(z^n),n = 0 .. infinity)

entsteht aus der L-Transformation durch Diskretisierung (Digitalisierung):

Tastet man stetige Funktionen f(t) mit einer Samplingrate a ab, so entstehen z.B. Treppenfunktionen

> digi:=(f,a)-> f(floor(a*t)/a);

digi := proc (f, a) options operator, arrow; f(floo...

> plot([10*sin(t),10*digi(sin,5),digi(x->x^2,10)/5],t=-2..7,numpoints=500);

[Maple Plot]

> #plot(sin(x),x=0..10,adaptive=false,sample=[seq(i,i=0..10)]); #so nicht :-))

>

Auf diese Treppenfunktionen läßt sich nun die L-Transformation

> F(s)=convert(laplace(f(t),t,s),int);

F(s) = int(f(t)*exp(-t*s),t = 0 .. infinity)

folgendermaßen anwenden: die Zeit t wird diskretisiert zu n*T, die L-Transformation innerhalb T für f(n)=const durchgeführt und aufsummiert, also z.B. für T=1

>

> sum(Int(f(n)*exp(-s*t),t=n..n+1),n=0..infinity)=sum(int(f(n)*exp(-s*t),t=n..n+1),n=0..infinity);

sum(Int(f(n)*exp(-t*s),t = n .. n+1),n = 0 .. infin...

oder (mit der Maus :-) vereinfacht:

> L:=(1-exp(-s))/s*sum(f(n)*exp(-s*n),n = 0 .. infinity);

L := (-exp(-s)+1)*sum(f(n)*exp(-s*n),n = 0 .. infin...

Dabei nennt man

> DL:=sum(f(n)*exp(-s*n),n = 0 .. infinity);

DL := sum(f(n)*exp(-s*n),n = 0 .. infinity)

die diskrete L-Transformation. Setzt man hierin exp(s)=z, so wird daraus

> subs(exp(-s*n)=1/z^n,DL);

sum(f(n)/(z^n),n = 0 .. infinity)

und zwischen der Z-Transformation und der L-Transformation besteht die Korrespondenz

> subs(exp(-s)=1/z,L/DL)*subs(exp(-s*n)=1/z^n,DL)*s=s*convert(laplace(f(t),t,s),int);

(-1/z+1)*sum(f(n)/(z^n),n = 0 .. infinity) = s*int(...

oder

> (1-1/z)*F(z)=s*F(s);

(-1/z+1)*F(z) = s*F(s)

Korespondenz heißt: es muss eine Treppenfunktion sein:

> (1-1/z)*ztrans(n,n,z)=s*inttrans[laplace](floor(t),t,s);

(-1/z+1)*z/((z-1)^2) = -exp(-s)/(-1+exp(-s))

> simplify(%);

1/(z-1) = -exp(-s)/(-1+exp(-s))

> simplify(subs(exp(-s)=1/z,%));

1/(z-1) = 1/(z-1)

und nicht eine stetige Funktion:

> simplify((1-1/z)*ztrans(n,n,z)=s*inttrans[laplace](t,t,s));

1/(z-1) = 1/s

Bei der Verwendung von ztrans ist zu beachten, dass z und F(z) im allgemeinen Fall komplex sind, und dass die Transformation nur definiert ist, wenn die Reihe konvergiert. Einfachstes Beispiel:

> ztrans(a,n,z);

a*z/(z-1)

Das ist der Grenzwert der geometrischen Reihe

> Sum(a/z^n,n=0..infinity)=sum(a/z^n,n=0..infinity);

Sum(a/(z^n),n = 0 .. infinity) = a*z/(-1+z)

der nur für |1/z|<1 existiert.

>

Die Rücktransformation gibt es natürlich auch:

> invztrans(rhs(%),z,n);

a

Weitere Beispiele

> ztrans( (-1)^n, n, z);

-z/(-z-1)

> simplify(ztrans(a^n,n,z));

-z/(-z+a)

> invztrans(%,z,n);

a^n

Einen Sonderfall (bzw. eine Analogie zur Diracfunktion) stellt die 'charakteristische Funktion' dar, die man erhält, wenn man eine Konstante im Bildbereich zurücktransformiert:

> invztrans(a,z,n);

charfcn[0](n)*a

Dabei meint charfcn[0](n) eine Folge, die nur für n=0 1 ist und Null sonst:

> seq(charfcn[0](n),n=-3..3);

0, 0, 0, 1, 0, 0, 0

und es gilt z.B.

> invztrans(1/z^5,z,n);

charfcn[5](n)

> ztrans(f(n)*charfcn[5](n),n,z);

f(5)/(z^5)

d.h., es wird aus der unendlichen Reihe der 5. Summand isoliert.

(Die charakteristische Funktion kann auch zum Morsen verwendet werden:

> seq(charfcn[-1,1,2](n),n=-3..4);

0, 0, 1, 0, 1, 1, 0, 0

> seq(charfcn[n](n),n=-3..4);

1, 1, 1, 1, 1, 1, 1, 1

> seq(charfcn[n mod 2](2*n mod 4),n=-3..20);

0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, ...

aber das nur nebenbei)

>

In Analogie zur L-Transformation besitzt die Z-Transformation die Eigenschaft, dass die Transformierte der Folge der Differenzen im wesentlichen das (z-1)-fache der Transformierten der Folge ist:

Erste Differenzen

> d1:=n->f(n+1)-f(n);

d1 := proc (n) options operator, arrow; f(n+1)-f(n)...

> alias(F(z)=ztrans(f(n),n,z)):

> collect(ztrans(d1(n),n,z),F(z));

(z-1)*F(z)-f(0)*z

Auf dieser Eigenschaft bauen die Lösungsverfahren für Differenzengleichungen auf

Kleiner Exkurs zu höheren Differenzen

Zweite Differenz

> d2:=n->d1(n+1)-d1(n);

d2 := proc (n) options operator, arrow; d1(n+1)-d1(...

> d2(k);

f(k+2)-2*f(k+1)+f(k)

> collect(ztrans(d2(n),n,z),F(z));

(-2*z+z^2+1)*F(z)+2*f(0)*z-f(0)*z^2-f(1)*z

Man kann das automatisieren:

a)

> d[1]:=n->f(n+1)-f(n):

> for i from 2 to 5
do d[i]:=unapply(d[i-1](n+1)-d[i-1](n),n); od:

> d[5](k);

f(k+5)-5*f(k+4)+10*f(k+3)-10*f(k+2)+5*f(k+1)-f(k)

b)

> dn:=n->add((-1)^(i+n)*binomial(n,i)*f(k+i),i=0..n);

dn := proc (n) options operator, arrow; add((-1)^(i...

> dn(5);

f(k+5)-5*f(k+4)+10*f(k+3)-10*f(k+2)+5*f(k+1)-f(k)

c)

> differ:=proc(n,k)
if n<2 then f(k+n)-f(k+n-1);
else differ(n-1,k+1)-differ(n-1,k) end if
end proc;

differ := proc (n, k) if n < 2 then f(k+n)-f(k+n-1)...
differ := proc (n, k) if n < 2 then f(k+n)-f(k+n-1)...

>

> differ(5,k);

f(k+5)-5*f(k+4)+10*f(k+3)-10*f(k+2)+5*f(k+1)-f(k)

> differ(5,1);

f(6)-5*f(5)+10*f(4)-10*f(3)+5*f(2)-f(1)

> seq(differ(i,k),i=1..5);

f(k+1)-f(k), f(k+2)-2*f(k+1)+f(k), f(k+3)-3*f(k+2)+...
f(k+1)-f(k), f(k+2)-2*f(k+1)+f(k), f(k+3)-3*f(k+2)+...
f(k+1)-f(k), f(k+2)-2*f(k+1)+f(k), f(k+3)-3*f(k+2)+...

Wenn z.B. die zweite Differenz konstant ist,

> ztrans(differ(2,n)=a,n,z);

>

z^2*F(z)-f(0)*z^2-f(1)*z-2*z*F(z)+2*f(0)*z+F(z) = a...

> solve(%,F(z));

z*(f(0)*z^2-3*f(0)*z+f(1)*z-f(1)+2*f(0)+a)/(z^3-3*z...

> collect(invztrans(%,z,n),n);

1/2*a*n^2+(f(1)-1/2*a-f(0))*n+f(0)

so muss die Folge quadratisch wachsen.

Geht man von den Differenzen zu den Differenzenquotienten über (Division durch die Schrittweite h^n), so kann man Differentialgleichungen durch Differenzengleichungen annähern, wichtigste Anwendung: finite Elemente.

>

>

Lösung von Differenzengleichungen mit der Z-Transformation

Erstes Beispiel

Die 'Schwingungsgleichung'

> alias(F(z)=ztrans(f(k),k,z)):

> dgl:=f(k+2)+f(k)=0;

dgl := f(k+2)+f(k) = 0

wird Z-transformiert

> bild:=ztrans(dgl,k,z);

bild := z^2*F(z)-f(0)*z^2-f(1)*z+F(z) = 0

> bildsol:=solve(bild,F(z));

bildsol := z*(f(0)*z+f(1))/(z^2+1)

nach der Transformierten aufgelöst und zurücktransformiert

> sol:=invztrans(bildsol,z,k);

sol := sum(1/2*((1/_alpha)^k*_alpha*f(0)-(1/_alpha)...

Lassen Sie sich nicht abschrecken - das funktioniert!

> sol:=eval(sol,{f(0)=0,f(1)=5});

sol := sum(-5/2*(1/_alpha)^k/_alpha,_alpha = RootOf...

Man kann es sogar zeichnen:

> plot(sol,k);

[Maple Plot]

Was haben wir da gezeichnet?

> add(allvalues(sol)[i],i=1..2);

5/2*I*(-I)^k-5/2*I*I^k

> evalc(Re(%));

5*sin(1/2*k*Pi)

Wir werden das in einem zweiten Beispiel näher untersuchen, zunächst der Vergleich mit rsolve:

> evalc(Re(rsolve(dgl,f)));

f(0)*cos(1/2*k*Pi)+f(1)*sin(1/2*k*Pi)

>

Zweites Beispiel

Paratstellen der Digitalisierungsfunktion und höherer Differenzen

> restart:

> digi:=(f,a)-> f(floor(a*k)/a):

> differ:=proc(n,k)
if n<2 then u(k+n)-u(k+n-1);
else differ(n-1,k+1)-differ(n-1,k) end if
end proc:

> alias(Delta=charfcn[0], U(z)=ztrans(u(k),k,z)):

>

Gedämpfte Schwingung

Die erste und zweite Differenz stehen nun als Differenzenquotienten mit der Schrittweite 1 für die erste und zweite Ableitung. Es können verschiedene Anregungen ('rechte Seiten') untersucht werden. Der schwarze Input läßt sich mit der rechten Maustaste zu 'Mapleinput' toggeln.

1. Ungestörtes System

> dgl:=8*differ(2,k)+4*differ(1,k)+3*u(k)=0;

dgl := 8*u(k+2)-12*u(k+1)+7*u(k) = 0

2. Einmalige Auslenkung

dgl:=8*differ(2,k)+4*differ(1,k)+3*u(k)=100;

3. Impuls zur Zeit Null

dgl:=8*differ(2,k)+4*differ(1,k)+3*u(k)=100*charfcn[0](k);

4. Sinusförmige Anregung

dgl:=8*differ(2,k)+4*differ(1,k)+3*u(k)=100*sin(2*Pi*k/3);

Transformation der Differenzengleichung und Auflösen nach der Bildfunktion

> bild:=ztrans( dgl, k, z);

bild := 8*z^2*U(z)-8*u(0)*z^2-8*u(1)*z-12*z*U(z)+12...

> Uz:=solve(bild, U(z));

Uz := 4*z*(2*u(0)*z+2*u(1)-3*u(0))/(8*z^2-12*z+7)

Rücktransformation

> uk:=invztrans( Uz, z, k);

uk := sum(1/5*(-2*(1/_alpha)^k*_alpha*u(0)+6*(1/_al...
uk := sum(1/5*(-2*(1/_alpha)^k*_alpha*u(0)+6*(1/_al...

>

Wer sich für _alpha interessiert, kann hier nachsehen

allvalues( op(selectfun(uk,RootOf)));

Zum Vergleich mit rsolve

rsolve(dgl,u);

>

Anfangswerte

> uka:=subs( u(0)=100, u(1)=50, uk);

uka := sum(1/5*(100*(1/_alpha)^k*_alpha+200*(1/_alp...

Kontrollausgabe

> seq(uka, k=0..10); #k:='k':

100, 50, -25/2, -125/2, -1325/16, -2225/32, -4075/1...

Plot

> plot(uka, k=0..20);

[Maple Plot]

Nanu - haben wir nun eine stetige Funktion?

> plot(uka, k=0..10,adaptive=false,sample=[seq(i,i=0..10)]);

[Maple Plot]

> ukf:=unapply(uka,k);

ukf := proc (k) options operator, arrow; sum(1/5*(1...

> plot([ukf(k),digi(ukf,1)],k=0..20);

[Maple Plot]

Also doch digital!

>

Neue Anfangswerte, neue Differenzengleichung

>

Zurück zum Inhaltsverzeichnis