Zurück zum Inhaltsverzeichnis

Kapitel 24

Reihenentwicklungen

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)];

typen := proc (x) options operator, arrow; [type(x,...

Taylor

> series(sin(x),x,10); typen(%);

series(1*x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^...

[true, true, true, series]

Laurent

> series(z/((z-(1+I))^2*(z-2)),z=1+I,2); typen(%);

series(-I*(z-1-I)^(-2)-I*(z-1-I)^(-1)+(1/2-1/2*I)+1...

[false, true, true, series]

> numapprox[laurent](z/((z-(1+I))^2*(z-2)),z=1+I,2); typen(%);

series(-I*(z-1-I)^(-2)-I*(z-1-I)^(-1)+(1/2-1/2*I)+1...

[false, true, true, series]

Taylor mit anderem Entwicklungspunkt

> taylor(z/((z-(1+I))^2*(z-2)),z=0,4); typen(%);

series(1/4*I*z+(1/4+3/8*I)*z^2+(1/2+3/16*I)*z^3+O(z...

[true, true, true, series]

Eine Laurentreihe...

> series(1/x,x=0);typen(%);

series(1*x^(-1),x)

[false, true, true, series]

Series, auch wenn ln(x) nicht etwickelt werden kann (um x=0)

> series(log(x),x); typen(%);

series(ln(x),x)

[false, false, true, series]

> series(log(x)+sin(x),x); typen(%);

series(ln(x)+1*x-1/6*x^3+1/120*x^5+O(x^7),x,7)

[false, false, true, series]

Entwicklungspunkt

> series(log(x)+sin(x),x=1,2); typen(%);

series(sin(1)+(1+cos(1))*(x-1)+O((x-1)^2),x=-(-1),2...

[true, true, true, series]

Kann nicht entwickelt werden (um x=0)

> series(sqrt(x),x=0); typen(%);

sqrt(x)

[false, false, false, `^`]

Verallgemeinerte Reihe

> series(sin(sqrt(x)),x,3); typen(%);

sqrt(x)-1/6*x^(3/2)+1/120*x^(5/2)+O(x^3)

[false, false, false, `+`]

> taylor(log(x),x);

Error, does not have a taylor expansion, try series()

> taylor(log(x),x=1,4);

series(1*(x-1)-1/2*(x-1)^2+1/3*(x-1)^3+O((x-1)^4),x...

>

Entwicklungspunkt infinity

> series(arctan(x), x=infinity);

1/2*Pi-1/x+1/3/(x^3)-1/5*1/(x^5)+O(1/(x^6))

Oder

> asympt(arctan(x),x);

1/2*Pi-1/x+1/3/(x^3)-1/5*1/(x^5)+O(1/(x^6))

> asympt(z/((z-(1+I))^2*(z-2)),z);

1/(z^2)+(4+2*I)/(z^3)+(8+10*I)/(z^4)+(8+28*I)/(z^5)...

> series(z/((z-(1+I))^2*(z-2)),z=infinity);

1/(z^2)+(4+2*I)/(z^3)+(8+10*I)/(z^4)+(8+28*I)/(z^5)...

>

Weiterverarbeitung von Reihen

Umwandlung in ein Polynom

> f:=series(sin(x), x=Pi);

f := series(-1*(x-Pi)+1/6*(x-Pi)^3-1/120*(x-Pi)^5+O...

> evalf(subs(x=2, f));

.9097898165+O((2-Pi)^6)

> whattype(f);

series

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);

p := -x+Pi+1/6*(x-Pi)^3-1/120*(x-Pi)^5

> evalf(subs(x=2, p));

.9097898165

> plot({p, sin(x)}, x=0..2*Pi);

[Maple Plot]

>

Es kann vorkommen, daß die Reihe abbricht

> f:=series(x^2, x=3);

f := series(9+6*(x-3)+1*(x-3)^2,x=-(-3))

trotzdem ist das Ergebnis vom Typ series

> whattype(f);

series

oder (C&P con f)

> f := series(9+6*(x-3)+1*(x-3)^2,x=-(-3));

f := series(9+6*(x-3)+1*(x-3)^2,x=-(-3))

und kann zunächst nicht vereinfacht werden

> simplify(f);

series(9+6*(x-3)+1*(x-3)^2,x=-(-3))

Wenn die 'exakte Reihe' zusammen mit anderen Ausdrücken vorkommt, wird sie eingeklammert

> f-x^2;

(series(9+6*(x-3)+1*(x-3)^2,x=-(-3)))-x^2

> simplify(convert(f,polynom)-x^2);

0

>

Rechnen mit Reihen

> s1:=series(sin(x),x=Pi/2);

s1 := series(1-1/2*(x-1/2*Pi)^2+1/24*(x-1/2*Pi)^4+O...

> s2:=series(cos(x), x=Pi/2);

s2 := series(-1*(x-1/2*Pi)+1/6*(x-1/2*Pi)^3-1/120*(...

> series(1/s1, x=Pi/2); # Reihe zu 1/sin(x)

series(1+1/2*(x-1/2*Pi)^2+5/24*(x-1/2*Pi)^4+O((x-1/...

> series(s1+x^2, x=Pi/2); # Reihe zu sin(x)+x^2

series((1+1/4*Pi^2)+Pi*(x-1/2*Pi)+1/2*(x-1/2*Pi)^2+...

und nicht

> s1+x^2;

(series(1-1/2*(x-1/2*Pi)^2+1/24*(x-1/2*Pi)^4+O((x-1...

>

> series(s1*s2, x=Pi/2); # Reihe zu sin(x)*cos(x)

series(-1*(x-1/2*Pi)+2/3*(x-1/2*Pi)^3-2/15*(x-1/2*P...

>

> diff(s1,x);

series(-1*(x-1/2*Pi)+1/6*(x-1/2*Pi)^3+O((x-1/2*Pi)^...

> int(s1,x);

series(1*(x-1/2*Pi)-1/6*(x-1/2*Pi)^3+1/120*(x-1/2*P...

Umkehrfunktion

> s1:=series(sin(x),x=0);

> series(arcsin(s1), x=0);

> simplify(convert(%, polynom));

s1 := series(1*x-1/6*x^3+1/120*x^5+O(x^6),x,6)

series(1*x+O(x^6),x,6)

x

>

>

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));

de := y(x)+diff(y(x),x) = sin(sqrt(x))

y(x) = exp(-x)*Int(sin(sqrt(u))*exp(u),u = 0 .. x)

Es gibt keine Stammfunktion

> integral:=int(sin(sqrt(x))*exp(x),x);

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);

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);

sr := sqrt(x)-1/6*x^(3/2)+1/120*x^(5/2)-1/5040*x^(7...
sr := sqrt(x)-1/6*x^(3/2)+1/120*x^(5/2)-1/5040*x^(7...

> de:=y(x)+diff(y(x),x)=sr:

> sol:=rhs(dsolve({de,y(0)=0},y(x)));

sol := 89909153/9909043200*x^(5/2)-750512033/396361...
sol := 89909153/9909043200*x^(5/2)-750512033/396361...
sol := 89909153/9909043200*x^(5/2)-750512033/396361...

Die Fehkerfunktion erf stört vielleicht etwas, also können wir nochmals in eine Reihe entwickeln

> psol:=convert(series(sol,x,9),polynom);

psol := 2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/...
psol := 2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/...

>

> plot([sr,sol,psol],x=0..100,-1..1,color=[red,blue,black]);

[Maple Plot]

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]);

[Maple Plot]

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]);

[Maple Plot]

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);

[Maple Plot]

> 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]);

[Maple Plot]

Also das System schwingt gedämpft mit. Allerdings scheint es Probleme zu geben

> plot([sr,psol],x=30..36,color=[red,blue]);

[Maple Plot]

Chaos durch Rundungsfehler?

> Digits:=100;

Digits := 100

> plot([sr,psol],x=30..38,color=[red,blue]);

[Maple Plot]

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]);

[Maple Plot]

> 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]);

[Maple Plot]

>

>

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));

de := y(x)+diff(y(x),x) = sin(sqrt(x))

y(x) = exp(-x)*Int(sin(sqrt(u))*exp(u),u = 0 .. x)

> integral:=int(sin(sqrt(x))*exp(x),x);

integral := int(sin(sqrt(x))*exp(x),x)

> integr:=diff(integral,x);

integr := sin(sqrt(x))*exp(x)

und integrieren

> inte:=int(series(integr,x,8),x);

inte := 2/3*x^(3/2)+1/3*x^(5/2)+41/420*x^(7/2)+461/...
inte := 2/3*x^(3/2)+1/3*x^(5/2)+41/420*x^(7/2)+461/...

mit der Exponentialfunktion multiplizieren

> se:=series(exp(-x),x,9);

se := series(1-1*x+1/2*x^2-1/6*x^3+1/24*x^4-1/120*x...

und das Ganze wieder in eine Reihe entwickeln. Lösung:

> los:=convert(series(se*inte,x,10),polynom);

los := 2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/2...
los := 2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/2...

Vielleicht sollte man die Probe machen?

> subs(y(x)=los,sin(sqrt(x))=sr,de);

2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/22680*x^...
2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/22680*x^...
2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/22680*x^...
2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/22680*x^...
2/3*x^(3/2)-1/3*x^(5/2)+41/420*x^(7/2)-493/22680*x^...

> eval(lhs(%)-rhs(%));

-3392923553/355687428096000*x^(17/2)

>

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);

siq := sqrt(x)-1/6*x^(3/2)+1/120*x^(5/2)-1/5040*x^(...
siq := sqrt(x)-1/6*x^(3/2)+1/120*x^(5/2)-1/5040*x^(...

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);

fy := a[1]*x^(3/2)+a[2]*x^(5/2)+a[3]*x^(7/2)+a[4]*x...

> fyd:=diff(fy, x);

fyd := 3/2*a[1]*sqrt(x)+5/2*a[2]*x^(3/2)+7/2*a[3]*x...

Die DGL lautet nun

> dgl:=fy+fyd-siq=0;

dgl := a[1]*x^(3/2)+a[2]*x^(5/2)+a[3]*x^(7/2)+a[4]*...
dgl := a[1]*x^(3/2)+a[2]*x^(5/2)+a[3]*x^(7/2)+a[4]*...
dgl := a[1]*x^(3/2)+a[2]*x^(5/2)+a[3]*x^(7/2)+a[4]*...

>

Und es sind die Koeffizienten a[i] zu bestimmen.

Etwas ordnen:

> collect(dgl,x);

(a[7]+1/1307674368000)*x^(15/2)+(15/2*a[7]+a[6]-1/6...
(a[7]+1/1307674368000)*x^(15/2)+(15/2*a[7]+a[6]-1/6...
(a[7]+1/1307674368000)*x^(15/2)+(15/2*a[7]+a[6]-1/6...

Zweite Reihe für Koeffizientenvergleich

> fb:=sum( b[i]*x^(i+1/2), i=0..7);

fb := b[0]*sqrt(x)+b[1]*x^(3/2)+b[2]*x^(5/2)+b[3]*x...

> match(lhs(dgl) = fb, x, 's');

true

> s;

{b[0] = -1+3/2*a[1], b[4] = -1/362880+a[4]+11/2*a[5...
{b[0] = -1+3/2*a[1], b[4] = -1/362880+a[4]+11/2*a[5...

>

Überzählige Gleichung entfernen und nur die rechten Seiten nehmen

> eqn:=map(rhs,remove(has,s,b[7]));

eqn := {15/2*a[7]+a[6]-1/6227020800, a[5]+1/3991680...
eqn := {15/2*a[7]+a[6]-1/6227020800, a[5]+1/3991680...

> sole:=solve(eqn);

sole := {a[7] = 757349/9340531200, a[2] = -1/3, a[1...

>

> assign(sole):

Vergleich von linker und rechter Seite der DGL

> plot({fy+fyd,siq},x=0..5,0..6);

[Maple Plot]

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);

coeftest := (a[7]+1/1307674368000)*x^(15/2)+(15/2*a...
coeftest := (a[7]+1/1307674368000)*x^(15/2)+(15/2*a...
coeftest := (a[7]+1/1307674368000)*x^(15/2)+(15/2*a...

Der nur für Polynome gedachte Befehl coeffs läßt sich nicht anwenden, aber

> factor(coeftest);

1/1307674368000*sqrt(x)*(1307674368000*x^7*a[7]+x^7...
1/1307674368000*sqrt(x)*(1307674368000*x^7*a[7]+x^7...
1/1307674368000*sqrt(x)*(1307674368000*x^7*a[7]+x^7...
1/1307674368000*sqrt(x)*(1307674368000*x^7*a[7]+x^7...
1/1307674368000*sqrt(x)*(1307674368000*x^7*a[7]+x^7...

>

> subs(sqrt(x)=1,x^7=0,factor(coeftest));

-1+15/2*x^6*a[7]+x^6*a[6]-1/6227020800*x^6+x^5*a[5]...
-1+15/2*x^6*a[7]+x^6*a[6]-1/6227020800*x^6+x^5*a[5]...

> eq:=coeffs(lhs(subs(x^7=0,subs(sqrt(x)=1,factor(coeftest)))),x);

eq := -1+3/2*a[1], 1/5040+9/2*a[4]+a[3], a[5]+1/399...
eq := -1+3/2*a[1], 1/5040+9/2*a[4]+a[3], a[5]+1/399...

> solve({%});

{a[7] = 757349/9340531200, a[2] = -1/3, a[1] = 2/3,...

>

4. Übrigens...

> de:=y(x)+diff(y(x),x)=sin(sqrt(x)):

> sollaplace:=dsolve(de,y(x),method=laplace);

sollaplace := y(x) = y(0)*exp(-x)+1/2*sqrt(Pi)*int(...
sollaplace := y(x) = y(0)*exp(-x)+1/2*sqrt(Pi)*int(...

> series(rhs(%),x);

Error, (in series/int) not implemented yet

So kommt man nicht weiter.

> #subs(sollaplace,de);

> #value(%);

>

Multivariable Taylorreihe

Beispiele

> mtaylor(log(x+y), [x=1, y=Pi], 3);

ln(1+Pi)+(x-1)/(1+Pi)+(y-Pi)/(1+Pi)-1/2*(x-1)^2/((1...

>

> f := sin(a+x)*cos(b+y);

f := sin(a+x)*cos(b+y)

> f1:=poisson(f, [x,y], 3 );

f1 := 1/2*sin(a+b)-1/2*sin(-a+b)+(-1/4*sin(a+b)+1/4...
f1 := 1/2*sin(a+b)-1/2*sin(-a+b)+(-1/4*sin(a+b)+1/4...
f1 := 1/2*sin(a+b)-1/2*sin(-a+b)+(-1/4*sin(a+b)+1/4...

> f2:=mtaylor(f, [x,y], 3);

f2 := sin(a)*cos(b)-sin(a)*sin(b)*y+cos(a)*x*cos(b)...
f2 := sin(a)*cos(b)-sin(a)*sin(b)*y+cos(a)*x*cos(b)...

> combine(f1-f2,trig);

0

> whattype(f1);

`+`

>

>

Formale Reihen

> restart: with(powseries);

[compose, evalpow, inverse, multconst, multiply, ne...
[compose, evalpow, inverse, multconst, multiply, ne...

>

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);

1, 1, 1/2, 1/6, 1/24, 1/120, 1/720

> p(100);

1/9332621544394415268169923885626670049071596826438...
1/9332621544394415268169923885626670049071596826438...

hier steht das Bildungsgesetz

> p(_k);

1/_k!

remember table

> op(4,op(p));

TABLE([0 = 1, 1 = 1, 2 = 1/2, 3 = 1/6, 4 = 1/24, 5 ...
TABLE([0 = 1, 1 = 1, 2 = 1/2, 3 = 1/6, 4 = 1/24, 5 ...
TABLE([0 = 1, 1 = 1, 2 = 1/2, 3 = 1/6, 4 = 1/24, 5 ...
TABLE([0 = 1, 1 = 1, 2 = 1/2, 3 = 1/6, 4 = 1/24, 5 ...

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);

series(1+1*x+1/2*x^2+1/6*x^3+1/24*x^4+O(x^5),x,5)

> whattype(%);

series

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);

series(1+1/2*x+1/2*x^2+1/4*x^3+1/8*x^4+O(x^5),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);

1, 0, -1/6, 0, 1/120, 0, -1/5040, 0, 1/362880, 0

> tpsform(mytayl,x);

series(cos(a)*x+(-1/2*sin(a))*x^2+(-1/6*cos(a))*x^3...

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);

series(1*x-1/6*x^3+1/120*x^5+O(x^6),x,6)

Oder allgemein

> eval(tpsform(mytayl,x),x=x-a);

series(cos(a)*(x-a)+(-1/2*sin(a))*(x-a)^2+(-1/6*cos...
series(cos(a)*(x-a)+(-1/2*sin(a))*(x-a)^2+(-1/6*cos...

>

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);

series(1*x-1/6*x^3+1/120*x^5-1/5040*x^7+O(x^9),x,9)...

> powcreate(q(k)=1/k!*signum(cos(Pi/2*k)));

> tpsform(q,x,9); #series(cos(x),x);

series(1-1/2*x^2+1/24*x^4-1/720*x^6+1/40320*x^8+O(x...

>

Nun kann mit den Reihen gerechnet werden

Produkt

> r:=evalpow(p*q);

r := proc (powparm) local nn, t1; option `Copyright...

> op(4,op(r));

TABLE([_k = sum(_powser1(j)*_powser2(_k-j),j = 0 .....

> #interface(verboseproc=3):print(r);interface(verboseproc=1):

> tpsform(r,x,9);

series(1*x-2/3*x^3+2/15*x^5-4/315*x^7+O(x^9),x,9)

> series( sin(x)*cos(x), x, 9);

series(1*x-2/3*x^3+2/15*x^5-4/315*x^7+O(x^9),x,9)

Umkehrung

> r:=reversion( p);

r := proc (powparm) local nn, t1; option `Copyright...

> tpsform(r, x, 9);

series(1*x+1/6*x^3+3/40*x^5+5/112*x^7+O(x^9),x,9)

> series( arcsin(x), x, 9);

series(1*x+1/6*x^3+3/40*x^5+5/112*x^7+O(x^9),x,9)

Ableitung

> r:=powdiff(p);

r := proc (powparm) local nn, t1; option `Copyright...

> tpsform(r,x,9);

series(1-1/2*x^2+1/24*x^4-1/720*x^6+1/40320*x^8+O(x...

>

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 );

r := proc (powparm) local nn, t1; option `Copyright...

> #op(4,op(r)); hier sind die 100 Koeffizienten

> t:=evalpow(p+r);

t := proc (powparm) local nn, t1; option `Copyright...

> t(81);

417527403333151351098563350293660679577328543206131...
417527403333151351098563350293660679577328543206131...
417527403333151351098563350293660679577328543206131...

> tpsform(t,x,15);

series(2*x+1/6*x^3+17/120*x^5+271/5040*x^7+7937/362...

> series(tan(x)+sin(x),x,15);

series(2*x+1/6*x^3+17/120*x^5+271/5040*x^7+7937/362...

>

Numerische Berechnung von Näherungsfunktionen

> restart:

> with(numapprox); with(orthopoly):

[chebdeg, chebmult, chebpade, chebsort, chebyshev, ...
[chebdeg, chebmult, chebpade, chebsort, chebyshev, ...

Zum Umgang mit with(orthopoly):

> T:=5;

T := 5

> T(2,x);

5

> with(orthopoly): T(2,x);

Warning, the name T has been redefined

2*x^2-1

> T;

T

D.h., die Warnung kommt _nach_ dem Überschreiben von T.

>

Funktion wählen

> f:=x->tan(x^2);

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

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);

[Maple Plot]

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;

mini0 := -.100575883+(4.867599399+(-36.71069895+(10...
mini0 := -.100575883+(4.867599399+(-36.71069895+(10...

-.100575883+4.867599399*x-36.71069895*x^2+102.48907...

.1006673752

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;

mini25 := (-.10387538e-2+(.74220446e-1+.9426059351*...
mini25 := (-.10387538e-2+(.74220446e-1+.9426059351*...

-(-.1038753800e-2+.7422044600e-1*x+.9426059351*x^2)...
-(-.1038753800e-2+.7422044600e-1*x+.9426059351*x^2)...

.5507722888e-3

> plot([f(x),mini0,mini25],x=-1/2..3,-10..10,colors,numpoints=500);

[Maple Plot]

>

Pade-Approximation (liefert hier mit dem geringsten Aufwand das beste Ergebnis, Enwicklung um den Ursprung)

> padef:=pade(f(x),x=a,[6,17]);

padef := (-10/99*x^6+x^2)/(1-43/99*x^4+17/1485*x^8+...

>

Tschebyscheff-Pade-Approximation

> chebpadef1:=chebpade(f(x),x=-x0+0.1..x0-0.1,[6,17]);

chebpadef1 := (.2485709044e-6+1.321535774*x^2-.2110...
chebpadef1 := (.2485709044e-6+1.321535774*x^2-.2110...
chebpadef1 := (.2485709044e-6+1.321535774*x^2-.2110...
chebpadef1 := (.2485709044e-6+1.321535774*x^2-.2110...

> #chebpadef1:=eval(%,T=orthopoly[T]);

>

> plot([f(x),padef,chebpadef1],x=-3..3,-10..10,colors,discont=true,numpoints=500);

[Maple Plot]

>

Zur Ergänzung noch die Entwicklung in Tschebyscheff-Reihen

> evalf(x0);

1.253314137

Vorsicht, das oszilliert!

> chebpadef2:=chebpade(f(x),x=a..1.25,[7,0]);

chebpadef2 := 6.035735130*T(0,1.600000000*x-1.00000...
chebpadef2 := 6.035735130*T(0,1.600000000*x-1.00000...
chebpadef2 := 6.035735130*T(0,1.600000000*x-1.00000...
chebpadef2 := 6.035735130*T(0,1.600000000*x-1.00000...
chebpadef2 := 6.035735130*T(0,1.600000000*x-1.00000...

> #chebpadef2;

> #chebpadef2:=eval(%,T=orthopoly[T]);

> cheby:=chebyshev(f(x),x=a..1,0.001);

cheby := .4863391011*T(0,2*x-1)+.6954069022*T(1,2*x...
cheby := .4863391011*T(0,2*x-1)+.6954069022*T(1,2*x...
cheby := .4863391011*T(0,2*x-1)+.6954069022*T(1,2*x...
cheby := .4863391011*T(0,2*x-1)+.6954069022*T(1,2*x...

> #cheby;

> #cheby:=eval(%,T=orthopoly[T]);

>

> plot([f(x),chebpadef2,cheby],x=-1/2..2,-10..10,colors,numpoints=500);

[Maple Plot]

>

Zurück zum Inhaltsverzeichnis