Kapitel 24
Für die Entwicklung von Funktionen in Reihen und vor allem für das Weiterarbeiten mit diesen Reihen ist es wichtig, sich mit den Typen vertraut zu machen. Natürlich kommt es auch auf den Entwicklungspunkt an.
> restart:
> typen:=x->[type(x,taylor),type(x,laurent),type(x,series),whattype(x)];
Taylor
> series(sin(x),x,10); typen(%);
Laurent
> series(z/((z-(1+I))^2*(z-2)),z=1+I,2); typen(%);
> numapprox[laurent](z/((z-(1+I))^2*(z-2)),z=1+I,2); typen(%);
Taylor mit anderem Entwicklungspunkt
> taylor(z/((z-(1+I))^2*(z-2)),z=0,4); typen(%);
Eine Laurentreihe...
> series(1/x,x=0);typen(%);
Series, auch wenn ln(x) nicht etwickelt werden kann (um x=0)
> series(log(x),x); typen(%);
> series(log(x)+sin(x),x); typen(%);
Entwicklungspunkt
> series(log(x)+sin(x),x=1,2); typen(%);
Kann nicht entwickelt werden (um x=0)
> series(sqrt(x),x=0); typen(%);
Verallgemeinerte Reihe
> series(sin(sqrt(x)),x,3); typen(%);
> taylor(log(x),x);
Error, does not have a taylor expansion, try series()
> taylor(log(x),x=1,4);
>
Entwicklungspunkt infinity
> series(arctan(x), x=infinity);
Oder
> asympt(arctan(x),x);
> asympt(z/((z-(1+I))^2*(z-2)),z);
> series(z/((z-(1+I))^2*(z-2)),z=infinity);
>
Weiterverarbeitung von Reihen
Umwandlung in ein Polynom
> f:=series(sin(x), x=Pi);
> evalf(subs(x=2, f));
> whattype(f);
Die Reihe muss erst in ein Polynom umgewandelt werden. (Alternative: Die Ausgabe von f mit C&P übernehmen und von Hand bearbeiten.)
> p:=convert(f, polynom);
> evalf(subs(x=2, p));
> plot({p, sin(x)}, x=0..2*Pi);
>
Es kann vorkommen, daß die Reihe abbricht
> f:=series(x^2, x=3);
trotzdem ist das Ergebnis vom Typ series
> whattype(f);
oder (C&P con f)
> f := series(9+6*(x-3)+1*(x-3)^2,x=-(-3));
und kann zunächst nicht vereinfacht werden
> simplify(f);
Wenn die 'exakte Reihe' zusammen mit anderen Ausdrücken vorkommt, wird sie eingeklammert
> f-x^2;
> simplify(convert(f,polynom)-x^2);
>
Rechnen mit Reihen
> s1:=series(sin(x),x=Pi/2);
> s2:=series(cos(x), x=Pi/2);
> series(1/s1, x=Pi/2); # Reihe zu 1/sin(x)
> series(s1+x^2, x=Pi/2); # Reihe zu sin(x)+x^2
und nicht
> s1+x^2;
>
> series(s1*s2, x=Pi/2); # Reihe zu sin(x)*cos(x)
>
> diff(s1,x);
> int(s1,x);
Umkehrfunktion
> s1:=series(sin(x),x=0);
> series(arcsin(s1), x=0);
> simplify(convert(%, polynom));
>
>
Differentialgleichungen mit Reihenentwicklungen lösen
Wir wollen testen, wie man mit Maple eine etwas ungewöhnliche DGL lösen kann.
> restart:
> de:=y(x)+diff(y(x),x)=sin(sqrt(x));
> dsolve({de,y(0)=0},y(x));
Es gibt keine Stammfunktion
> integral:=int(sin(sqrt(x))*exp(x),x);
>
und deshalb auch keine Reihe
> series(integral,x);
Error, (in series/int) unable to compute series
> integralu:=int(sin(sqrt(u))*exp(u),u=0..x);
> series(integralu,x);
Error, (in series/int) not implemented yet
>
Nun gibt es mehrere Möglichkeiten
1. Die rechte Seite (Inhomogenität) in eine Reihe entwickeln
> sr:=convert(series(sin(sqrt(x)),x,9),polynom);
> de:=y(x)+diff(y(x),x)=sr:
> sol:=rhs(dsolve({de,y(0)=0},y(x)));
Die Fehkerfunktion erf stört vielleicht etwas, also können wir nochmals in eine Reihe entwickeln
> psol:=convert(series(sol,x,9),polynom);
>
> plot([sr,sol,psol],x=0..100,-1..1,color=[red,blue,black]);
Eine höhere Ordnung für die Entwicklung der Fehlerfunktion
> psol:=convert(series(sol,x,100),polynom):
> plot([sr,sol,psol],x=0..200,-1..1,color=[red,blue,black]);
Wenn man also nicht gezwungen ist, die Fehlerfunktion in eine Reihe zu entwickeln, sollte man es lassen!
> plot([sr,sol,sin(sqrt(x))],x=0..200,-2..2,color=[red,blue,black]);
Aber auch die ursprüngliche Entwicklung von sin(sqrt(x)) ist verbesserungsfähig - hier wird ja das CAS erst interessant!
> restart:
> sr:=convert(series(sin(sqrt(x)),x,30),polynom):
> plot([sin(sqrt(x)),sr],x=0..600,color=[red,black],numpoints=500);
> de:=y(x)+diff(y(x),x)=sr:
Die DGL mit dem Reihenansatz für die rechte Seite lösen
> sol:=rhs(dsolve( {de,y(0)=0}, y(x))):
und wieder in eine Reihe entwickeln (in diesem Fall geht series/int). Interessant dabei ist, dass man wegen exp zu ziemlich hohen Ordnungen greifen muß (aber bitte nicht übertreiben und das Worksheet vorher abspeichern).
> psol:=convert(series(sol,x,100),polynom):
> plot([sr,psol],x=0..36,color=[red,blue]);
Also das System schwingt gedämpft mit. Allerdings scheint es Probleme zu geben
> plot([sr,psol],x=30..36,color=[red,blue]);
Chaos durch Rundungsfehler?
> Digits:=100;
> plot([sr,psol],x=30..38,color=[red,blue]);
Es ist eben schwierig, eine 'Gerade' durch ein Polynom anzunähern - und das auch noch bei x=38. Aber ansonsten kann sich die Lösung innerhalb einer Schwingung schon sehen lassen.
> plot([sr,psol],x=0..38,color=[red,blue]);
> Digits:=10:
>
Die Lösung mit der nicht entwickelten Fehlerfunktion ist dagegen (je nach Rechengenauigkeit) auch in einem 10-fach größeren Bereich stabil.
> assume(x>0):
> resol:=evalc(Re(value(sol))):
> plot([sr,resol],x=0..450,color=[red,blue]);
>
>
2. Den Integranden in der von Maple vorgeschlagenen Lösung in eine Reihe entwickeln.
> restart:
> sr:=convert(series(sin(sqrt(x)),x,9),polynom):
> de:=y(x)+diff(y(x),x)=sin(sqrt(x));
> dsolve({de,y(0)=0},y(x));
> integral:=int(sin(sqrt(x))*exp(x),x);
> integr:=diff(integral,x);
und integrieren
> inte:=int(series(integr,x,8),x);
mit der Exponentialfunktion multiplizieren
> se:=series(exp(-x),x,9);
und das Ganze wieder in eine Reihe entwickeln. Lösung:
> los:=convert(series(se*inte,x,10),polynom);
Vielleicht sollte man die Probe machen?
> subs(y(x)=los,sin(sqrt(x))=sr,de);
> eval(lhs(%)-rhs(%));
>
3. Schließlich kann man es mit einer benutzerdefinierten Reihe für die Lösung selbst versuchen (wenn alle Stricke reissen...).
Zunächst wieder die Reihe für die rechte Seite
> siq:=convert(series( sin(sqrt(x)), x, 8),polynom);
Und nun die Reihen für die linke Seite mit unbekannten Koeffizienten
> a:=array(1..7):
> fy:=sum( a[i]*x^(i+1/2), i=1..7);
> fyd:=diff(fy, x);
Die DGL lautet nun
> dgl:=fy+fyd-siq=0;
>
Und es sind die Koeffizienten a[i] zu bestimmen.
Etwas ordnen:
> collect(dgl,x);
Zweite Reihe für Koeffizientenvergleich
> fb:=sum( b[i]*x^(i+1/2), i=0..7);
> match(lhs(dgl) = fb, x, 's');
> s;
>
Überzählige Gleichung entfernen und nur die rechten Seiten nehmen
> eqn:=map(rhs,remove(has,s,b[7]));
> sole:=solve(eqn);
>
> assign(sole):
Vergleich von linker und rechter Seite der DGL
> plot({fy+fyd,siq},x=0..5,0..6);
Nachdem wir nur in niedriger Ordnung gerechnet haben, war eigentlich auch nicht mehr zu erwarten.
>
>
Noch eine Variante für den Zugriff auf die Koeffizienten
> a:='a':
> coeftest:=collect(dgl,x);
Der nur für Polynome gedachte Befehl coeffs läßt sich nicht anwenden, aber
> factor(coeftest);
>
> subs(sqrt(x)=1,x^7=0,factor(coeftest));
> eq:=coeffs(lhs(subs(x^7=0,subs(sqrt(x)=1,factor(coeftest)))),x);
> solve({%});
>
4. Übrigens...
> de:=y(x)+diff(y(x),x)=sin(sqrt(x)):
> sollaplace:=dsolve(de,y(x),method=laplace);
> series(rhs(%),x);
Error, (in series/int) not implemented yet
So kommt man nicht weiter.
> #subs(sollaplace,de);
> #value(%);
>
Beispiele
> mtaylor(log(x+y), [x=1, y=Pi], 3);
>
> f := sin(a+x)*cos(b+y);
> f1:=poisson(f, [x,y], 3 );
> f2:=mtaylor(f, [x,y], 3);
> combine(f1-f2,trig);
> whattype(f1);
>
>
> restart: with(powseries);
>
Mit powcreate wird eine Prozedur erzeugt, die die Koeffizienten berechnet und in der remember table ablegt
> powcreate( p(n)=1/n!);
>
> seq(p(i),i=0..6);
> p(100);
hier steht das Bildungsgesetz
> p(_k);
remember table
> op(4,op(p));
Ausgabe der Prozedur, falls erwünscht
> #interface(verboseproc=3):print(p);interface(verboseproc=1):
>
Die so gebildeten Koeffizienten können in der 'truncated powerseries form' dargestellt werden
> tpsform(p, x, 5);
> whattype(%);
Bevor eine neue Reihe mit dem alten Namen gebildet werden soll, sollte man den Namen wieder freigeben
> p:='p':
Die Koeffizienten können auch rekursiv und mit Anfangsbedingungen angegeben werden
> powcreate( p(n)=p(n-1)*p(n-2), p(0)=1, p(1)=1/2);
> tpsform(p, x, 5);
>
Grundsätzlich kann man sich nun seine eigene Taylorreihe aufstellen
> mytayl:='mytayl':
>
powcreate(mytayl(n)=diff(sin(a),a$(n))/(n)!,mytayl(0)=0):
> seq(eval(mytayl(i),a=0),i=1..10);
> tpsform(mytayl,x);
Das geht nicht:
> tpsform(mytayl,x-a);
Error, (in tpsform) wrong number (or type) of parameters in function series
> eval(tpsform(mytayl,x),a=0);
Oder allgemein
> eval(tpsform(mytayl,x),x=x-a);
>
Wenn man die Taylorformel nicht verwendet, muss man unter Umständen dafür sorgen, dass Koeffizienten Null werden
> p:='p':
> powcreate(p(k)=1/k!*signum(sin(Pi/2*k)));
> tpsform(p,x,9); #series(sin(x),x);
> powcreate(q(k)=1/k!*signum(cos(Pi/2*k)));
> tpsform(q,x,9); #series(cos(x),x);
>
Nun kann mit den Reihen gerechnet werden
Produkt
> r:=evalpow(p*q);
> op(4,op(r));
> #interface(verboseproc=3):print(r);interface(verboseproc=1):
> tpsform(r,x,9);
> series( sin(x)*cos(x), x, 9);
Umkehrung
> r:=reversion( p);
>
tpsform(r, x, 9);
> series( arcsin(x), x, 9);
Ableitung
> r:=powdiff(p);
> tpsform(r,x,9);
>
Die Koeffizienten kann man sich auch anders besorgen: Die Übernahme von Reihen, die mit series gebildet wurden, erfolgt mit powpoly. Dabei kann man leicht 'beliebig viele' Koeffizienten zur Verfügung stellen
> s:=series( tan(x), x=0, 100):
> #op(s); # hier ist die Struktur von series
> r:=powpoly( convert(s,polynom), x );
> #op(4,op(r)); hier sind die 100 Koeffizienten
> t:=evalpow(p+r);
> t(81);
> tpsform(t,x,15);
> series(tan(x)+sin(x),x,15);
>
Numerische Berechnung von Näherungsfunktionen
> restart:
> with(numapprox); with(orthopoly):
Zum Umgang mit with(orthopoly):
> T:=5;
> T(2,x);
> with(orthopoly): T(2,x);
Warning, the name T has been redefined
> T;
D.h., die Warnung kommt _nach_ dem Überschreiben von T.
>
Funktion wählen
> f:=x->tan(x^2);
Unstetigkeitsstelle
> x0:=sqrt(Pi/2):
Bereichsangaben
> a:=0: b:=1:
Zunächst eine Reihenentwicklung um x0 (erzeugt eine Laurentreihe)
> sx0:=series(f(x),x=x0,40):
Taylorreihe um den Ursprung
> s:=series(f(x),x=a,40):
> colors:=color=[red,blue,black]:
> plot([f(x),convert(s,polynom),convert(sx0,polynom)],x=-2..3,-10..10,colors,numpoints=500,legend=["tan(x^2)","Reihe um 0","Reihe um x0"],discont=true);
Die Laurentreihe ist am Pol nicht von der Funktion zu unterscheiden, hat aber für kleine x eine sehr große Abweichung von f(x). Die Taylorreihe approxiemiert die Funktion für kleine x sehr gut.
>
Minimax approximiert nach dem Remes-Algorithmus. Wenn man bei der ersten Variante (Polynom) mit der rechten Bereichsgrenze zu nahe an x0 geht, beginnt die Näherung zu oszillieren.
> mini0:=minimax(f(x),x=a..x0-0.1,5,1,'fehler');expand(mini0);fehler;
Mit einer gebrochen rationalen Funktion erhält man eine wesentlich bessere Näherung
> mini25:=minimax(f(x),x=a..x0-0.05,[2,5],1,'maxerror');normal(expand(mini25));maxerror;
> plot([f(x),mini0,mini25],x=-1/2..3,-10..10,colors,numpoints=500);
>
Pade-Approximation (liefert hier mit dem geringsten Aufwand das beste Ergebnis, Enwicklung um den Ursprung)
> padef:=pade(f(x),x=a,[6,17]);
>
Tschebyscheff-Pade-Approximation
> chebpadef1:=chebpade(f(x),x=-x0+0.1..x0-0.1,[6,17]);
> #chebpadef1:=eval(%,T=orthopoly[T]);
>
> plot([f(x),padef,chebpadef1],x=-3..3,-10..10,colors,discont=true,numpoints=500);
>
Zur Ergänzung noch die Entwicklung in Tschebyscheff-Reihen
> evalf(x0);
Vorsicht, das oszilliert!
> chebpadef2:=chebpade(f(x),x=a..1.25,[7,0]);
> #chebpadef2;
> #chebpadef2:=eval(%,T=orthopoly[T]);
> cheby:=chebyshev(f(x),x=a..1,0.001);
> #cheby;
> #cheby:=eval(%,T=orthopoly[T]);
>
> plot([f(x),chebpadef2,cheby],x=-1/2..2,-10..10,colors,numpoints=500);
>