Kapitel 18
Symbolisches Lösen von Differentialgleichungen
Ohne Anfangsbedingungen
> restart:
Formulierung mit dem D-Operator
>
dsolve(D(y)(x) + y(x) = 1, y(x));
oder mit diff
> dsolve( diff(y(x), x) + y(x) = 1, y(x));
Mit Anfangsbedingung für die Funktion
> de := diff(y(x),x) * y(x) * (1+x^2) = x;
> dsolve( { de, y(0)=0 }, y(x));
Probe
> subs(%[1], de);
> simplify(%);
>
Mit Anfangsbedingung für die Ableitung
> dsolve( { de, D(y)(1)=1 }, y(x));
>
Implizite Lösung
> dsolve( (y(x)^2 - x)*D(y)(x) + x^2-y(x) = 0, y(x) ,implicit);
Explizite Lösung (default)
> dsolve( (y(x)^2 - x)*D(y)(x) + x^2-y(x)=0, y(x), explicit );
>
>
Lösungsbasis und partikuläre Lösung
> de:=diff(y(x), x$4) + 5*diff(y(x),x$2) - 36*y(x)=sin(x);
>
> simplify(dsolve( de, y(x)));
>
>
> simplify(dsolve( de, y(x), output=basis));
>
>
>
Beispiel1: Freier Fall
> restart:
> sol:=dsolve({ m*D(D(z))(t) = -m*g, z(0)=10, D(z)(0)=0 }, z(t));
Übernahme der Lösung mit assign(). (Nur einmal möglich: restart.)
> assign(sol);
>
> z(t);
Es wird nur z(t) zugewiesen:
> z(x);
>
Neue Funktion
> funcz:=unapply(rhs(sol), t);
> funcz(x);
>
Werte
> g:=9.81;
> solve( funcz(t)=0, t);
> t0:=max(%);
> #plot( {funcz(t), diff(funcz(t),t)}, t=0..t0);
>
Beispiel 2: Federpendel mit Dämpfung
> restart;
> m:=1: c:=1/3: k:=1:
> dsolve( {m*D(D(y))(t) + c*D(y)(t) + k*y(t)=0, y(0)=1, D(y)(0)=0}, y(t));
> funcy:=rhs(%);
> plot(funcy, t=0..18);
>
Platzhalter für Lösungen: DESol
Diese DGL ist nicht geschlossen lösbar (Fadenpendel):
> restart:
> dgl:=diff(x(t),t$2) = -g/l*sin(x(t));
>
dsolve({dgl, x(0)=0, D(x)(0)=v0}, x(t));
Warning, computation interrupted
>
Reihenentwicklung
> #Order:=10:
> reihe:=dsolve({dgl, x(0)=x0, D(x)(0)=v0}, x(t),series);
>
> dgsol:=DESol( dgl, x(t), {x(0)=x0, D(x)(0)=v0} );
>
> intdg:=int(dgsol,t);
> diff(intdg,t,t,t)+g/l*sin(diff(intdg,t));
>
> intreihe:=int(rhs(reihe),t);
> diff(intreihe,t,t,t)+g/l*sin(diff(intreihe,t));
> fehler:=convert(%,polynom);
> collect(%,x0);
>
> g,l,x0,v0:='g','l','x0','v0';
> g,l,x0,v0:=1,1,0.1,0;
> fehler;
> plot(fehler,t=-10..10);
>
>
Systeme von Differentialgleichungen
Ein typisches System wie es z.B. beim Radioaktiven Zerfall vorkommt, zunächst allgemein formuliert:
> restart:
> sys:=D(x)(t)=a11*x(t)+a12*y(t),D(y)(t)=a21*x(t)+a22*y(t);
>
> dsolve( {sys}, {x(t), y(t)} );
>
Mit spezieller Wahl der Koeffizienten und Anfangsbedingungen
> a11,a12,a21,a22:=-1,0,1,-0.4;
> sol:=dsolve( {sys, x(0)=1, y(0)=0}, {x(t), y(t)} );
> funcx:=subs(sol, x(t));
> funcy:=subs(sol, y(t));
> plot( [funcx, funcy], t=0..6);
>
Oder im Phasendiagramm
> plot( [funcx, funcy, t=0..4]);
>
Lösung durch Laplace-Transformation
Sprungantwort
> restart:
> alias(sigma=Heaviside):
> u:=sigma(t-1)- sigma(t-5):
> de:=diff(y(t), t) + y(t) = u:
> dsolve( {de}, y(t));
> dsolve( {de}, y(t), method=laplace);
>
> sol1 := rhs(%);
> simplify( subs(y(t)=sol1, de) );
> lhs(%)-rhs(%);
> y(0):=0:
> plot( {u, sol1}, t=0..10);
>
Näherungslösung durch Reihenentwicklung
> restart:
> Order:=10:
> de:=diff(y(x), x$2) + diff(y(x),x) + y(x) = x+sin(x);
> sol1:=dsolve( {de, y(0)=0, D(y)(0)=0}, y(x), series);
> f1:=convert( rhs(sol1), polynom);
> sol2:=dsolve( {de, y(0)=0, D(y)(0)=0}, y(x));
> f2:=rhs(sol2):
> plot( {f1,f2}, x=0..2*Pi, 0..8);
>
Mit infolevel die Lösungsmethoden anzeigen:
> restart:
> infolevel[dsolve]:=3:
> dsolve( 3*D(y)(x) + y(x)* (1 + (2*x-1)*y(x)^3) = 0, y(x));
Methods for first order ODEs:
Trying to isolate the derivative dy/dx...
Successful isolation of dy/dx
-> Trying classification methods
trying a quadrature
trying 1st order linear
trying Bernoulli
Bernoulli successful
>
>
>
Höhere Ordnung
> alias(y=y(x)):
> dsolve( diff(y, x$4) + 5*diff(y,x$2) - 36*y=0, y);
Methods for high order ODEs:
Trying to isolate the derivative d^4y/dx^4...
Successful isolation of d^4y/dx^4
-> Trying classification methods
trying a quadrature
trying high order exact linear fully integrable
trying linear constant coefficient
linear constant coefficient successful
> infolevel[dsolve]:=1:
>
Numerische Lösung von Differentialgleichungen
Beispiel 1
Zum Vergleich die wird zunächst die symbolische Lösung bereitgestellt (gedämpfte Schwingung):
> restart:
> de:=m*diff( y(t), t$2) +c*diff(y(t),t)+k*y(t)=0;
> sol1:=dsolve( {de, y(0)=1, D(y)(0)=0}, y(t) ):
> funcy1:=unapply(rhs(sol1),t):
> m:=1: c:=1: k:=1:
Numerische Auswertung
> evalf(funcy1(3));
>
Numerische Lösung
> funcy2:=dsolve( {de, y(0)=1, D(y)(0)=0}, y(t), numeric );
> funcy2(3);
> subs( funcy2(3), y(t));
Wertetabelle
> seq( [ t=n*0.2, y=subs( funcy2(n*0.2), y(t)) ], n=0..10 ):
> array( [%] );
> Digits;
> with(plots):
Warning, the name changecoords has been redefined
> odeplot( funcy2, [t,y(t)], 0..10);
> plot( funcy1(x)-'subs(funcy2(x),y(t))',x=0..10);
>
Beispiel 2
Proton im Magnetfeld
> restart: with(linalg): with(plots):
Warning, the protected names norm and trace have been redefined and unprotected
Warning, the name changecoords has been redefined
Daten
> m:= 1.6e-27: # Masse [kg]
> q:= 1.6e-19: # Ladung [Coulomb]
> bz:= 1: # Magnetfeld in z-Richtung [Tesla]
> r:= [x(t), y(t), z(t)]: # Ortsvektor
> v:= map(diff, r, t): # Geschwindigkeitsvektor
> a:= map(diff, v, t): # Beschleunigungsvektor
> b:= vector( [0,0,bz] ): # Feldvektor
> sys := evalm( m*a - q*crossprod(v,b));
> de:=seq( sys[i]=0, i=1..3);
>
Für glatte Kurven sollte die Rechengenauigkeit vor dem Aufruf von dsolve erhöht werden (insbesondere bei dem gewählten Beispiel mit hohen Zehnerpotenzen).
> Digits:=11:
> conditions:= x(0)=0, y(0)=0, z(0)=0, D(x)(0)=1e8, D(y)(0)=0, D(z)(0)=1e6:
Es gibt in Maple zwei Arten der Ausgabe der erzeugten Prozeduren:
procedurelist erzeugt eine Prozedur, die die Lösung als liste ausgibt (default):
> sol:=dsolve( {de, conditions}, convert(r,set), numeric);
Mit den auskommentierten Befehlen wird die Prozedur angezeigt.
> #interface(verboseproc=3):
> #print(sol); interface(verboseproc=1):
> solp:=dsolve( {de, conditions}, convert(r,set), numeric,output=procedurelist);
> solp(1e-8);
>
listprocedure erzeugt eine Liste von Prozeduren:
> soll:=dsolve( {de, conditions}, convert(r,set), numeric,output=listprocedure);
> soll(1e-8);
>
Die Liste der Prozeduren in den Ortsvektor einsetzen:
> rs:=subs(soll, r);
> rs(1e-8); # Test
>
Zeichnen
> spacecurve( [ 'rs[1](t)','rs[2](t)','rs[3](t)' ],
> t=0..3e-7,numpoints=500,axes=framed,color=black);
>
Die Komponenten des Ortsvektors müssen unausgewertet übergeben werden:
> eval(rs[1](t),t=1e-7);
Error, (in rs[1]) cannot evaluate boolean: abs(t)-abs(-.299398796030000024e-6+t) < 0
>
>
Grafische Darstellung von Differentialgleichungen
odeplot
Numerische Lösungen mit odeplot darstellen
> restart: with(plots):
Warning, the name changecoords has been redefined
> sol:=dsolve( { diff(y(t),t)+3*y(t)=sin(t)^2, y(0)=0 }, y(t), numeric):
>
Mit der Option numpoints bekommt man glatte Kurven.
> odeplot( sol, [t, y(t)], 0..10,numpoints=500);
>
System von Gifferentialgleichungen.
> sol:=dsolve(
> { diff(y(t),t)=x(t)+sin(t), diff(x(t),t)=-y(t), x(0)=0, y(0)=0 },
> { x(t), y(t) }, numeric);
>
> odeplot( sol, [t, x(t)], 0..20,numpoints=500);
> odeplot( sol, [t, y(t)], 0..20,numpoints=500);
>
>
2D-Phasendiagramm
> odeplot( sol, [x(t), y(t)], 0..20,
> view=[-10..10, -10..10],
> numpoints=500,
> scaling=constrained);
3D-Phasendiagramm
> odeplot( sol, [t, x(t), y(t)], 0..20, shading=none, axes=boxed, numpoints=500);
>
DEtools
Richtungsfelder mit dfieldplot
> restart: with(DEtools):
> dfieldplot( diff(y(t),t)=sin(y(t)), [y(t)], t=-1..5, y=-1..4);
>
dfieldplot( [diff(x(t),t)=sin(x(t))+y(t), diff(y(t),t)=sin(y(t))-x(t)],
[x(t),y(t)], t=0..1, x=-1..1, y=-1..1);
Richtungsfelder und Lösungskurven zu gegebenen Anfangsbedingungen mit phaseportrait:
>
phaseportrait( diff(y(t),t)=sin(y(t)), [y(t)],
t=-1..10, [[y(4)=1], [y(4)=2], [y(4)=-1], [y(4)=-2]],linecolor=black);
>
Ohne Richtungsfeld
> start:=seq( [0,0,n*0.1], n=0..10);
> phaseportrait( [diff(x(t),t)=sin(x(t))+y(t), diff(y(t),t)=sin(y(t))-x(t)],
> [x(t),y(t)], t=0..10, {start}, stepsize=0.1, arrows=NONE, linecolor=black);
>
Richtungsfeld mit DEplot
> alias(x=x(t), y=y(t)):
> DEplot( diff(x,t)=1+x-x^2, [x], t=0..3,
> {seq([0,n/2], n=-3..5)}, x=-1..3,
> dirgrid=[10,10], arrows=LARGE);
>
Gleichungssystem
> DEplot( [diff(x,t)=-y, diff(y,t)=x-y], [x,y], -1..1, { [0,0,1] },linecolor=black );
>
3D-Phasendiagramm
> DEplot3d( [diff(x,t)=sin(y)-x/10, diff(y,t)=x], [x,y], 0..40, {[0,1,1]},
> stepsize=0.1, linecolor=black );
> alias(x=x,y=y):
> de := diff(y(x), x$4) + diff(y(x), x) = sin(x);
>
Schar von Lösungskurven
>
DEplot(de, [y(x)], x=0..2,
[seq( [ y(1)=1, D(y)(1)=-1, (D@@2)(y)(1)=1, (D@@3)(y)(1)=n ], n=-5..5)],
linecolor=black );
>
Partielle Differentialgleichungen
Erster Ordnung
>
restart:
alias(f=f(x,y)):
> pde:= x*diff(f,x) + y*diff(f,y)=f;
> sol:= pdsolve( pde, f);
> subs(sol, pde);
> simplify(lhs(%)-rhs(%));
>
> pde:=diff(f,x)+diff(f,y)=2*(x+y)*f;
> sol:= pdsolve( pde, f);
> subs(sol, pde);
> simplify(lhs(%)-rhs(%));
>
Mit RootOf
> pde:=x^4*diff(f,x)^2 + y^2*diff(f,x)^2=f;
> sol:=pdsolve( pde, f);
> subs(sol, pde);
> simplify(lhs(%)-rhs(%));
>
Zweiter Ordnung
>
restart:
alias(u=u(t,x)):
Wellengleichung
> pde:=diff(u,t$2)=a^2*diff(u,x$2);
> sol:=pdsolve(pde,u);
> subs(sol,pde): simplify(lhs(%)-rhs(%));
> pde:=diff(u,t$2)=a^2*diff(u,x$2)+sin(t);
> sol:=pdsolve(pde,u);
> subs(sol,pde): simplify(lhs(%)-rhs(%));
>
Wärmeleitung :
> pde:=diff(u,x$2)-diff(u,t)=0;
Nur ein Vorschlag
> sol:=pdsolve(pde,u);
Vergleiche Worksheet zu Integraltransformationen.
>
>
restart:
alias(u=u(t,x,y)):
2D Wellengleichung
> pde:=diff(u,t$2) = a^2 * (diff(u,x$2)+diff(u,y$2));
> sol:=pdsolve(pde,u);
>
Laplace Gleichung : ok
> pde:=diff(u,x$2)+diff(u,y$2)=0;
> pdsolve(pde, u);
>
Höhere Ordnung
> restart:
> alias(f=f(x,y));
Keine Lösung
> #infolevel[all]:=0:
>
#infolevel[all]:=9: # Vorsicht: langer Ausdruck
pdsolve( diff(f,y) + f*diff(f,x)+diff(f,x$3)=0, f );
>
>
>
Partielle DGln, grafische Darstellung
Aus dem reichhaltigen Angebot von PDEplot ein paar Beispiele
> restart: with(PDEtools);
>
> pd:=2*diff(z(x,y),x)-z(x,y)*diff(z(x,y),y)=1/10;
Parameterdarstellung [cos(s),sin(s),s]
> PDEplot( pd , z(x,y),[cos(s),sin(s),s], [s=-6..6], axes=framed,lightmodel='light4',numchar=60,color=yellow,orientation=[50,80],x=-4..4,y=-4..4);
>
Mit Charakteristiken (basechar=true)
> PDEplot( pd, z(x,y),[cos(s),sin(s),s], [s=-12..12], axes=framed,lightmodel='light4',color=yellow,numchar=60,orientation=[8,60],numsteps=[-20,40], stepsize=0.05,basechar=true);
>
Animation
> PDEplot( pd, z(x,y),[cos(s),sin(s),s],[s=-12..12],axes=framed,lightmodel='light4',color=yellow,numchar=60,orientation=[40,85],animate=true);
>
>
>