kapitel-29maple.mw

 

 

Mathematik  

 

von T. Arens, F. Hettlich, Ch. Karpfinger, U. Kockelkorn, K. Lichtenegger,  H. Stachel  

 

(zu Kapitel 29:  Partielle Differenzialgleichungen) 

 

 

> restart;
 

>
 

Der Befehl pdsolve 

 

Auch zum Loesen partieller Differentialgleichungen stellt  Maple einen Befehl und das Programmpaket PDEtools  bereit. Zentral ist der allgemeine Befehl   pdsolve , der in speziellen Situationen explizit Loesungen finden kann.  Versuchen wir zunaechst eine quasi-lineare Differenzialgleichung erster Ordnung. 

> pde:= x^2*diff(u(x,y),x)+x*y*diff(u(x,y),y)-(u(x,y))^2=0;
pdsolve(pde);
 

 

`+`(`*`(`^`(x, 2), `*`(diff(u(x, y), x))), `*`(x, `*`(y, `*`(diff(u(x, y), y)))), `-`(`*`(`^`(u(x, y), 2)))) = 0
u(x, y) = `/`(`*`(x), `*`(`+`(1, `*`(_F1(`/`(`*`(y), `*`(x))), `*`(x)))))
 

Die vom System eingefuehrte Funktion _F1 ist dabei frei waehlbar. Mit den Optionen zum pdsolve Kommando bieten sich eine ganze Reihe von Moeglichkeiten  partielle Differenzialgleichung zu untersuchen, bis hin zu numerischen Loesungen (Option numeric)  durch die Methode der finiten Differenzen.  Darueberhinaus gibt es das spezielle Paket PDEtools, aus dem wir uns noch ein paar Befehle ansehen.   

> with(PDEtools);
 

[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ, ConservedCurrentTest, ConservedCurrents, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet, InfinitesimalGenerator, Infinitesimals, IntegratingFactorT...
 

Wir koennen die Ordnung der Differenzialgleichung erfragen 

> difforder(pde);
 

1
 

oder uns das charakteristische System verschaffen 

> charstrip(pde,u(x,y));
 

{diff(u(_s), _s) = `*`(`^`(u(_s), 2)), diff(x(_s), _s) = `*`(`^`(x(_s), 2)), diff(y(_s), _s) = `*`(x(_s), `*`(y(_s)))}
 

Wichtig ist die Moeglichkeit Variablentransformationen vorzunehmen, z.B. eine Transformation, wie wir es bei der d'Alembert Loesung zur Waermeleitungsgleichung im Abschnitt 29.2 gesehen haben. Dazu dient der Befehl dchange  

> dchange({x=t+s,y=t-s},pde);
 

`+`(`*`(`^`(`+`(t, s), 2), `*`(`+`(`*`(`/`(1, 2), `*`(diff(u(s, t), s))), `*`(`/`(1, 2), `*`(diff(u(s, t), t)))))), `*`(`+`(t, s), `*`(`+`(t, `-`(s)), `*`(`+`(`-`(`*`(`/`(1, 2), `*`(diff(u(s, t), s)))...
`+`(`*`(`^`(`+`(t, s), 2), `*`(`+`(`*`(`/`(1, 2), `*`(diff(u(s, t), s))), `*`(`/`(1, 2), `*`(diff(u(s, t), t)))))), `*`(`+`(t, s), `*`(`+`(t, `-`(s)), `*`(`+`(`-`(`*`(`/`(1, 2), `*`(diff(u(s, t), s)))...
 

Gibt es weitere Informationen zur Loesung u, so lassen sich diese durch dsubs auch in die Differenzialgleichung einsetzen. 

> eq :=u(x,y)=(v(x,y))^2;
dsubs(eq,pde);
 

 

u(x, y) = `*`(`^`(v(x, y), 2))
`+`(`*`(2, `*`(`^`(x, 2), `*`(v(x, y), `*`(diff(v(x, y), x))))), `*`(2, `*`(x, `*`(y, `*`(v(x, y), `*`(diff(v(x, y), y)))))), `-`(`*`(`^`(v(x, y), 4)))) = 0
 

Auf diesem Weg koennten wir auch eine Separation ansetzen. 

> eq :=u(x,y)=v(x)*w(y);
dsubs(eq,pde);
 

 

u(x, y) = `*`(v(x), `*`(w(y)))
`+`(`*`(`^`(x, 2), `*`(diff(v(x), x), `*`(w(y)))), `*`(x, `*`(y, `*`(v(x), `*`(diff(w(y), y))))), `-`(`*`(`^`(v(x), 2), `*`(`^`(w(y), 2))))) = 0
 

Ob eine additive oder multiplikative  Separation sinnvoll ist, laesst sich mit dem Befehl separability erfragen. So laesst sich die Laplacegleichung sowohl durch eine Summe als auch durch ein Produkt separieren  

> pde:=diff(u(x,y),x,x) + diff(u(x,y),y,y)=0;
separability(pde,u(x,y));
separability(pde,u(x,y),`*`);
separability( lhs(pde)+(u(x,y))^2=0,u(x,y));
 

 

 

 

`+`(diff(diff(u(x, y), x), x), diff(diff(u(x, y), y), y)) = 0
0
0
`+`(`*`(2, `*`(diff(u(x, y), y), `*`(diff(u(x, y), x)))))
 

Das Resultat 0 erhalten wir, wenn eine Separation vollstaendig moeglich ist. Dagegen ist das letzte Ergebnis wie folgt zu interpretieren: Fuer die modifizierte nichtlineare Differenzialgleichung kann nur unter der Bedingung, dass der angegebene Ausdruck verschwindet, eine Loesung in der separierten Form `+`(f(x), g(y)) gefunden werden kann. 

   

Uebrigens, wenn Ihnen die Indexschreibweise fuer die partiellen Ableitungen lieber ist, so bietet Maple mit dem Befehl difftable eine entsprechende Moeglichkeit. 

>
 

Charakteristiken  

> restart; with(plots):
 

Neben den oben beschriebenen allgemeinen Befehlen zum Loesen partieller Differenzialgleichungen, koennen wir selbstverstaendlich auch die einzelnen Schritte der im Buch vorgestellten Methoden mit Maple nachvollziehen. Wir betrachten in diesem Abschnitt das Charakteristikenverfahren und im naechsten einen Separationsansatz im Detail. 

 

Starten wir mit der quasi linearen partiellen Differenzialgleichung erster Ordnung (siehe Abschnitt 29.3). Zunaechst geben wir die Differenzialgleichung ein. 

> pde:= D[1](u)(x,y) - 2*u(x,y)*D[2](u)(x,y) - u(x,y) = 0;
 

`+`((D[1](u))(x, y), `-`(`*`(2, `*`(u(x, y), `*`((D[2](u))(x, y))))), `-`(u(x, y))) = 0
 

und lesen das charakteristische System ab, 

> ode1:= D(k1)(s) = 1; ode2:= D(k2)(s) = 2*w(s); ode3:= D(w)(s) = w(s);
 

 

 

(D(k1))(s) = 1
(D(k2))(s) = `+`(`*`(2, `*`(w(s))))
(D(w))(s) = w(s)
 

(oder beschaffen es uns mit charstrip).  Die allgemeine Loesung des Systems gewoehnlicher Differenzialgleichungen ergibt sich zu 

> lsgtmp:=dsolve([ode1,ode2,ode3]);
 

{k1(s) = `+`(s, _C3), k2(s) = `+`(`*`(2, `*`(_C2, `*`(exp(s)))), _C1), w(s) = `*`(_C2, `*`(exp(s)))}
 

 

Bemerkung: An dieser Stelle ergibt sich ein technisches Problem, dass Ihnen vielleicht nicht auffaellt. Aber die Zuordnungen der Konstanten und die Reihenfolge der Angaben der Funktionen in lsgtmp kann variieren, sodass wir im Worksheet einwenig Aufwand treiben muessen, um die richtigen Kombinationen zu extrahieren. Wenn Sie die Ausgabe vor Augen haben, koennen Sie k1, k2, w einfach durch eine Zeile wie  

   k1:= rhs(lsg_tmp[1]); k2:=rhs(lsg_tmp[2]); w:=rhs(lsg_tmp[3]); 

definieren. Wir muessen hier aber erst die Position der Eintraege in der Loesungsmenge festlegen und bestimmen. 

> lsgtmp:=convert(lsgtmp,list);
ltmp:=map(lhs,lsgtmp):
member('k1(s)',ltmp,'pk1'): member('k2(s)',ltmp,'pk2'): member('w(s)',ltmp,'pw'):
k1:=rhs(lsgtmp[pk1]);
k2:=rhs(lsgtmp[pk2]);
w:=rhs(lsgtmp[pw]);
 

 

 

 

[k1(s) = `+`(s, _C3), k2(s) = `+`(`*`(2, `*`(_C2, `*`(exp(s)))), _C1), w(s) = `*`(_C2, `*`(exp(s)))]
`+`(s, _C3)
`+`(`*`(2, `*`(_C2, `*`(exp(s)))), _C1)
`*`(_C2, `*`(exp(s)))
 

 

Nun stellen wir uns die  Anfangskurve durch gamma(t) = (0, t) parametrisiert vor und erhalten fuer die Konstanten _C1,_C2,_C3 in Abhaengigkeit von t an der ausgewaehlten Stelle s=0 die Bedingungen  

> bed1:= subs(s=0,k1)=0;
bed2:= subs(s=0,k2)=t;
bed3:=subs(s=0,w)=2*subs(s=0,k2);
 

 

 

_C3 = 0
`+`(`*`(2, `*`(_C2, `*`(exp(0)))), _C1) = t
`*`(_C2, `*`(exp(0))) = `+`(`*`(4, `*`(_C2, `*`(exp(0)))), `*`(2, `*`(_C1)))
 

und bekommen 

> lsgtmp:=solve({bed1,bed2,bed3},[_C1,_C2,_C3]);
 

[[_C1 = `+`(`-`(`*`(3, `*`(t)))), _C2 = `+`(`*`(2, `*`(t))), _C3 = 0]]
 

Wir setzen diese Werte ein, sodass sich  

> lsgtmp:=[seq(lsgtmp[1,j],j=1..3)];
ltmp:=map(lhs,lsgtmp):
member('_C1',ltmp,'c1'): member('_C2',ltmp,'c2'): member('_C3',ltmp,'c3'):
_C1:=rhs(lsgtmp[c1]);
_C2:=rhs(lsgtmp[c2]);
_C3:=rhs(lsgtmp[c3]);
 

 

 

 

[_C1 = `+`(`-`(`*`(3, `*`(t)))), _C2 = `+`(`*`(2, `*`(t))), _C3 = 0]
`+`(`-`(`*`(3, `*`(t))))
`+`(`*`(2, `*`(t)))
0
 

Wieder haben wir dabei  zunaechst die Reihenfolge der Variablen in der Loesung ermittelt. Nun loesen wir das Gleichungssystem  

> lsgtmp:=solve({x=k1,y=k2},[s,t]);
 

[[s = x, t = `/`(`*`(y), `*`(`+`(`-`(3), `*`(4, `*`(exp(x))))))]]
 

und setzen diese Ausdruecke in w ein um die Loesung u zu bestimmen 

> lsgtmp:=[seq(lsgtmp[1,j],j=1..2)];
ltmp:=map(lhs,lsgtmp):
member('s',ltmp,'ps'): member('t',ltmp,'pt'):
u:=subs(s=rhs(lsgtmp[ps]),t=rhs(lsgtmp[pt]),w);
 

 

[s = x, t = `/`(`*`(y), `*`(`+`(`-`(3), `*`(4, `*`(exp(x))))))]
`+`(`/`(`*`(2, `*`(y, `*`(exp(x)))), `*`(`+`(`-`(3), `*`(4, `*`(exp(x)))))))
 

Eine Illustration der Loesung, ermoeglicht uns Chrakteristiken und Grundcharakteristiken einzuzeichnen  

> t:=-1/4:
p1:=plot3d(u,x=-0.1..1,y=-1..1,axes=boxed,style=patchnogrid):
p2:=spacecurve([k1,k2,0],s=0..0.5,color=black,thickness=2):
p3:=spacecurve([k1,k2,2*exp(s)*t],s=-0.1..0.5,color=red,thickness=2):
p4:=spacecurve([0,y,0],y=-1..1,color=blue,thickness=2):
t:=0:
p5:=spacecurve([0,y,2*y],y=-1..1,color=yellow,thickness=2):
p6:=spacecurve([k1,k2,2*t*exp(s)],s=-0.1..1,color=red,thickness=2):
t:=1/4:
p7:=spacecurve([k1,k2,0],s=-0.1..0.5,color=black,thickness=2):
p8:=spacecurve([k1,k2,2*t*exp(s)],s=-0.1..0.5,color=red,thickness=2):
display(p1,p3,p4,p5,p6,p7,p8);
 

Plot
 

Eingezeichnet in den Graphen der Loesung u sind einige  Charakteristiken (rot), eine zugehoerige Grundcharakteristik (schwarz), die Kurve, auf der Anfangsbedingungen gegeben sind (blau) und die Kurve mit den Anfangswerten (gelb). Bei den 3D Bildern hilft es oft, das Bild einwenig zu drehen, um es sich besser vortellen zu koennen. 

>
 

Ein Paukenschlag 

 

> restart: with(plots): with(DEtools):
 

Wir nutzen die Gelegenheit um noch ein paar weitere Befehle aus dem DEtool Paket kennenzulernen.  Als Beispiel einer Separation versuchen wir die Schwingung einer kreisfoermigen Membran zu analysieren. Die zugehoerige Schwingungsgleichung lautet 

> pde:= D[1,1](u)(t,x,y) - D[2,2](u)(t,x,y) - D[3,3](u)(t,x,y) = 0;
 

`+`((D[1, 1](u))(t, x, y), `-`((D[2, 2](u))(t, x, y)), `-`((D[3, 3](u))(t, x, y))) = 0
 

wobei wir in diesem Beispiel den Differentialoperator D nutzen, also im Gegensatz zu oben u als Funktion und nicht als Ausdruck behandeln. Nehmen wir an, das die Membran fest eingespannt ist, so ergibt sich eine Randbedingung von der Form  

 u(t,x,y)  = 0     fuer    `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = 1.     

 

Eine Separation der Zeitabhaengigkeit bedeutet 

    u(t, x, y) = `*`(w(t), `*`(v(x, y))) 

Um die partielle Differentialgleichung in diesen Variablen zu schreiben nutzen wir den Befehl Dchangevar im Paket DEtools und erhalten 

> pde:=Dchangevar(u(t,x,y)=w(t)*v(x,y),pde,t);
 

`+`(`*`(((`@@`(D, 2))(w))(t), `*`(v(x, y))), `-`(`*`(w(t), `*`((D[1, 1](v))(x, y)))), `-`(`*`(w(t), `*`((D[2, 2](v))(x, y))))) = 0
 

Wir teilen durch w und v, um die Abhaengigkeiten zu separieren, 

> pde:=collect(pde/w(t)/v(x,y),D);
 

`+`(`-`(`/`(`*`((D[1, 1](v))(x, y)), `*`(v(x, y)))), `-`(`/`(`*`((D[2, 2](v))(x, y)), `*`(v(x, y)))), `/`(`*`(((`@@`(D, 2))(w))(t)), `*`(w(t)))) = 0
 

und erhalten  zwei DGL'n, die  durch  eine Konstante gekoppelt sind: 

> dgl_t:=remove(has,op(1,pde),x)=-k^2;
dgl_xy:=remove(has,op(1,pde),t)=k^2;
 

 

`/`(`*`(((`@@`(D, 2))(w))(t)), `*`(w(t))) = `+`(`-`(`*`(`^`(k, 2))))
`+`(`-`(`/`(`*`((D[1, 1](v))(x, y)), `*`(v(x, y)))), `-`(`/`(`*`((D[2, 2](v))(x, y)), `*`(v(x, y))))) = `*`(`^`(k, 2))
 

 

Fuer die Zeitabhaengigkeit w ergibt sich eine gewoehnliche Differentialgleichung zweiter Ordnung, deren allgemeine Loesung, eine harmonische Schwingung,  wir leicht ermitteln koennen  

> lsg_t:=dsolve({dgl_t,w(0)=0},w(t));
 

w(t) = `*`(_C1, `*`(sin(`*`(k, `*`(t)))))
 

wobei wir hier die Anfangsbedingung w(0) = 0 verwenden.  

 

Nun bleibt noch die Ortsabhaengigkeit, also die  Helmholtzgleichung 

> dgl_xy:= dgl_xy*v(x,y);
 

`*`(v(x, y), `*`(`+`(`-`(`/`(`*`((D[1, 1](v))(x, y)), `*`(v(x, y)))), `-`(`/`(`*`((D[2, 2](v))(x, y)), `*`(v(x, y))))))) = `*`(v(x, y), `*`(`^`(k, 2)))
 

zu untersuchen.  Da wir Loesungen im Einheitskreis suchen waehlen wir Polarkoordinaten und transformieren die partielle Differntialgleichung mit PDEchangecoords auf die neuen Koordinaten.  

> dgl_p:=PDEchangecoords(dgl_xy,[x,y],polar,[r,phi]);
 

`/`(`*`(`+`(`-`(`*`((D[1](v))(r, phi), `*`(r))), `-`((D[2, 2](v))(r, phi)), `-`(`*`((D[1, 1](v))(r, phi), `*`(`^`(r, 2)))), `-`(`*`(v(`*`(r, `*`(cos(phi))), `*`(r, `*`(sin(phi)))), `*`(`^`(k, 2), `*`(...
 

> dgl_p:=subs(v(r*cos(phi),r*sin(phi))=v(r,phi),dgl_p);
 

`/`(`*`(`+`(`-`(`*`((D[1](v))(r, phi), `*`(r))), `-`((D[2, 2](v))(r, phi)), `-`(`*`((D[1, 1](v))(r, phi), `*`(`^`(r, 2)))), `-`(`*`(v(r, phi), `*`(`^`(k, 2), `*`(`^`(r, 2))))))), `*`(`^`(r, 2))) = 0
 

 

Nun folgt eine Separationsansatz in den Polarkoordinaten, d.h.                     

    v(r, phi) = rad(r) arg(phi), 

> dgl_p:=Dchangevar(v(r,phi)=rad(r)*arg(phi),dgl_p,r);
 

`/`(`*`(`+`(`-`(`*`((D(rad))(r), `*`(arg(phi), `*`(r)))), `-`(`*`(rad(r), `*`(((`@@`(D, 2))(arg))(phi)))), `-`(`*`(((`@@`(D, 2))(rad))(r), `*`(arg(phi), `*`(`^`(r, 2))))), `-`(`*`(rad(r), `*`(arg(phi)...
`/`(`*`(`+`(`-`(`*`((D(rad))(r), `*`(arg(phi), `*`(r)))), `-`(`*`(rad(r), `*`(((`@@`(D, 2))(arg))(phi)))), `-`(`*`(((`@@`(D, 2))(rad))(r), `*`(arg(phi), `*`(`^`(r, 2))))), `-`(`*`(rad(r), `*`(arg(phi)...
`/`(`*`(`+`(`-`(`*`((D(rad))(r), `*`(arg(phi), `*`(r)))), `-`(`*`(rad(r), `*`(((`@@`(D, 2))(arg))(phi)))), `-`(`*`(((`@@`(D, 2))(rad))(r), `*`(arg(phi), `*`(`^`(r, 2))))), `-`(`*`(rad(r), `*`(arg(phi)...
 

Wir sortieren den Ausdruck  

> dgl_p:=collect(r^2*dgl_p/rad(r)/arg(phi),D);
 

`+`(`-`(`/`(`*`(r, `*`((D(rad))(r))), `*`(rad(r)))), `-`(`/`(`*`(((`@@`(D, 2))(arg))(phi)), `*`(arg(phi)))), `-`(`/`(`*`(`^`(r, 2), `*`(((`@@`(D, 2))(rad))(r))), `*`(rad(r)))), `-`(`*`(`^`(k, 2), `*`(...
 

und extrahieren die beiden separierten Differntialgleichungen. Zunaechst bezueglich phi  

> dgl_phi:=remove(has,op(1,dgl_p),r)=c;
 

`+`(`-`(`/`(`*`(((`@@`(D, 2))(arg))(phi)), `*`(arg(phi))))) = c
 

Es ergibt sich die allgemeine Loesung 

> lsg_phi:=dsolve(dgl_phi,arg(phi));
 

arg(phi) = `+`(`*`(_C1, `*`(sin(`*`(`^`(c, `/`(1, 2)), `*`(phi))))), `*`(_C2, `*`(cos(`*`(`^`(c, `/`(1, 2)), `*`(phi))))))
 

> arg:=unapply(rhs(lsg_phi),phi);
 

proc (phi) options operator, arrow; `+`(`*`(_C1, `*`(sin(`*`(`^`(c, `/`(1, 2)), `*`(phi))))), `*`(_C2, `*`(cos(`*`(`^`(c, `/`(1, 2)), `*`(phi)))))) end proc
 

Wir benoetigen alle  `+`(`*`(2, `*`(Pi))) -periodischen Loesungen, die wir durch 

> _EnvAllSolutions:=true:
lsg_c:=solve({arg(0)=arg(2*Pi),D(arg)(0)=D(arg)(2*Pi)},c);
 

{c = 0}, {c = `*`(`^`(_Z1, 2))}
 

bekommen, wobei an die vom Sytem eingefuehrte Hilfsvariable _Z die Bedingung  

> about(indets(subs(lsg_c[2],c)));
 

{_Z1}:
 

 is used in the following assumed objects
 [_Z1] assumed integer
 

zu stellen ist.  Wir fuehren statt c den Index n ein, um diese Loesungen zu kennzeichnen 

> assume(n,integer); assume(n>=0); c:=n^2;
 

`*`(`^`(n, 2))
 

Also erhalten wir die Loesungen 

> lsg_phi:=simplify(arg(phi));
 

`+`(`*`(_C1, `*`(sin(`*`(n, `*`(phi))))), `*`(_C2, `*`(cos(`*`(n, `*`(phi))))))
 

 

Als letztes fehlt uns noch die radiale Abhaengigkeit. Wir betrachten also die Differentialgleichung fuer die Funktion rad(r) 

> dgl_r:=simplify(eval(dgl_p))*rad(r);
 

`+`(`-`(`*`(r, `*`((D(rad))(r)))), `*`(`^`(n, 2), `*`(rad(r))), `-`(`*`(`^`(r, 2), `*`(((`@@`(D, 2))(rad))(r)))), `-`(`*`(`^`(k, 2), `*`(`^`(r, 2), `*`(rad(r)))))) = 0
 

> dgl_r:=collect(dgl_r,rad(r));
 

`+`(`*`(`+`(`*`(`^`(n, 2)), `-`(`*`(`^`(k, 2), `*`(`^`(r, 2))))), `*`(rad(r))), `-`(`*`(r, `*`((D(rad))(r)))), `-`(`*`(`^`(r, 2), `*`(((`@@`(D, 2))(rad))(r))))) = 0
 

Diese gewoehnliche Differentialgleichung heisst Besselsche Differentialgleichung und als allgemeine Loesung ergeben sich Linearkombinationen von Bessel- und Neumannfunktion n-ter Ordnung  

> lsg_rad:=dsolve(dgl_r,rad(r));
 

rad(r) = `+`(`*`(_C1, `*`(BesselJ(n, `*`(k, `*`(r))))), `*`(_C2, `*`(BesselY(n, `*`(k, `*`(r))))))
 

 

Zur Illustration schauen wir uns einige Graphen dieser Funktionen an, deren Auswertung durch die Befehle BesselJ und BesselY  implementiert ist. 

> J0:= plot(BesselJ(0,t), t=0..10, color=red):
J1:= plot(BesselJ(1,t), t=0..10, color=blue):
J2:= plot(BesselJ(2,t), t=0..10, color=green):
display(J0,J1,J2);
 

Plot_2d
 

und die Neumannfunktionen 

> Y0:= plot(BesselY(0,t), t=0..10, y=-1..1, color=red):
Y1:= plot(BesselY(1,t), t=0..10, y=-1..1, color=blue):
Y2:= plot(BesselY(2,t), t=0..10, y=-1..1, color=green):
display(Y0,Y1,Y2);
 

Plot_2d
 

 

Wir sehen, dass die Neumannfunktionen bei r = 0 singulaer werden. Dies kann uns keine sinnvolle Beschreibung der Schwingung der Membran liefern, sodass wir diese Loesungen ausschliessen koennen. 

> assume(k>0):
lsg_rad:= rhs(subs(_C2=0,lsg_rad));
 

`*`(_C1, `*`(BesselJ(n, `*`(k, `*`(r)))))
 

 

Es fehlt nun noch die Randbedingung bei r = 1 zu beruecksichtigen.   Aus  u =0, d.h. lsg_rad(1) = 0, ergibt sich, dass fuer die Kopplungskonstante k nur Nullstellen der Besselfunktion in frage kommen. 

> solve(subs(r=1,lsg_rad)=0,k);
 

RootOf(BesselJ(n, _Z))
 

Diese Nullstellen sind in Maple durch das Kommando BesselJZeros fest verankert und wir koennen uns einige Werte beschaffen. 

> z01:=evalf(BesselJZeros(0,1));
z11:=evalf(BesselJZeros(1,1));
z21:=evalf(BesselJZeros(2,1));
z02:=evalf(BesselJZeros(0,2));
z12:=evalf(BesselJZeros(1,2));
 

 

 

 

 

2.404825558
3.831705970
5.135622302
5.520078110
7.015586670
 

 

Mit diesen Nullstellen kann die allgemeine Loesung durch eine Reihe 

    u(t, r, phi) = sum(`*`(`+`(`*`(a[n], `*`(sin(`*`(m, `*`(phi))))), `*`(b[n], `*`(cos(`*`(m, `*`(phi)))))), `*`(J(m, `*`(z[n], `*`(r))), `*`(sin(`*`(z[n], `*`(t)))))), n = 1 .. infinity) 

mit  den der Groesse nach nummerierten Nullstellen z[n] der zugehoerigen m-ten Besselfunktionen angegeben werden. Jeder Summand repraesentiert einen Schwingungsmodus (eine Eigenschwingung) und deren Ueberlagerung liefert den Klang, wobei die Anteile, die Koeffizienten, der einzelnen Modi von der Anregung der Schwingung abhaengen. Betrachten wir eine kleine Auswahl dieser Moden 

> u01:= BesselJ(0,z01*r)*sin(z01*t);
u21:= cos(2*phi)*BesselJ(2,z21*r)*sin(z21*t);
u12:= cos(phi)*BesselJ(1,z12*r)*sin(z12*t);
 

 

 

`*`(BesselJ(0, `+`(`*`(2.404825558, `*`(r)))), `*`(sin(`+`(`*`(2.404825558, `*`(t))))))
`*`(cos(`+`(`*`(2, `*`(phi)))), `*`(BesselJ(2, `+`(`*`(5.135622302, `*`(r)))), `*`(sin(`+`(`*`(5.135622302, `*`(t)))))))
`*`(cos(phi), `*`(BesselJ(1, `+`(`*`(7.015586670, `*`(r)))), `*`(sin(`+`(`*`(7.015586670, `*`(t)))))))
 

und lassen mit animate3d einen Modus schwingen  

> animate3d([r*cos(phi), r*sin(phi), u12], r=0..1, phi=0..2*Pi, t=1..3.7, frames=40);
 

Plot
 

oder eine Ueberlagerung der Schwingungen. 

> animate3d([r*cos(phi), r*sin(phi), 0.05*u01+0.05*u21+0.1*u12], r=0..1, phi=0..2*Pi, t=1..3.7, frames=50);
 

Plot
 

>
 

Aufgaben 

1. Loesen Sie Aufgabe 29.7. mit Maple, d.h. gesucht ist die Loesung des Anfangsrandwertproblems  

     `+`(diff(u(x, t), x, x), `*`(4, `*`(diff(u(x, t), t))), `-`(`*`(3, `*`(u(x, t))))) = 0 

mit den Bedingungen u(x, 0) = `*`(`^`(sin(x), 3)) und u(0, t) = 0,

 

Loesung 

>
 

. 

2. Bestimmen Sie die Loesung u des Anfangswertproblems 

  `+`(`*`(x, `*`(diff(u(x, y), x))), `/`(`*`(x, `*`(diff(u(x, y), y))), `*`(y, `*`(u(x, y)))), u(x, y)) = 0 

und u(`*`(`^`(t, 2)), t) = 1 fuer `<`(0, t) (s. Aufagbe 29.11). 

Loesung 

>