Differentiaalvergelijkingen in Mathematica

Hier illustreren we hoe we in Mathematica met differentiaalvergelijkingen kunnen werken. We doen dit aan de hand van het voorbeeld 
y'  =  (E^(-x/3))y+1, met y(0)=-2.

Het lijnelementenveld

We kunnen het lijnelementenveld van een differentiaalvergelijking
tekenen met behulp van de functie PlotVectorField. Deze bevindt
zich in het package PlotField, dat dus eerst geladen moet worden:

<<Graphics`PlotField`

(Zie het notebook `laden.nb' over wat er mis kan gaan als dit commando
niet op tijd uitgevoerd wordt.)
Laat  y(x)  een oplossing van de differentiaalvergelijking zijn. Dan heeft
de raaklijn aan y(x) in het punt (x,y(x)) als richtings-vector (1, y'(x)).
Nu kennen we  y(x) niet, maar  y'(x) wel (uit de differentiaalvergelijking).
We kunnen dus het lijnelementenveld van de differentiaalvergelijking
tekenen door in punten (x,y) de vector (1, y' ) te tekeken (waar y' met
behulp van de differentiaalvergelijking uitgerekend wordt).   
Het commando PlotVectorField doet dat. In het onderstaande voorbeeld
zijn alleen de eerste drie argumenten van wezenlijk belang: het eerste argument
geeft de vector die getekend moet worden, en het tweede en derde argument
respectievelijk de x-range, en y-range. De andere argumenten worden 
beneden nog kort toegelicht. Om het commando te gebruiken hoef je het
gebruik van deze argumenten echter niet te begrijpen; je kunt ze eenvoudigweg
overnemen.

plot1 = PlotVectorField[ { 1, E^(-x/3) y +1 }, {x,0,20}, {y,-10,10 },
Axes -> True, ScaleFunction -> (1.0&), HeadWidth -> 0 ]

Hier zijn de opties Axes -> True (om de assen te krijgen), ScaleFunction -> (1.0&) (om alle lijntjes even lang te krijgen), en HeadWidth -> 0 
(om te voorkomen dat Mathematica pijltjes tekent in plaats van lijntjes)
gebruikt. Ga na wat er gebeurt als je een of meerdere van deze opties
weglaat.

De vergelijking oplossen

We kunnen een vergelijking oplossen met de functie DSolve. Deze functie
kan de algemene oplossing van een differentiaalvergelijking vinden, of 
de oplossing met een gegeven beginvoorwaarde.  De algemene oplossing
vinden we door aan DSolve drie argumenten te geven: de vergelijking, 
de onbekende functie, en tenslotte de onafhankelijke variabele. In het 
benedenstaande voorbeeld is  y[x]  de onbekende functie, en x de onafhankelijke 
variabele.

sol = DSolve[ y'[x] == E^(-x/3) y[x] + 1, y[x], x ]

Hier merken we drie dingen op: in de oplossing staat een functie die we niet 
kennen, ExpIntegralEi. Mathematica kent veel meer standaardfuncties dan wij.
Op dit ogenblik is het niet zo van belang wat de functie ExpIntegralEi precies
is; want bij het gebruik maakt dat niet uit. Ten tweede zien we dat er in de
oplossing een onbekende konstante voorkomt, namelijk C[1]. Dit komt
doordat we de algemene oplossing berekend hebben. Ten derde is de oplossing
gegeven als een substitutieregel; (of preciezer: als een lijst van substitutieregels,
hier ter lengte 1), als we de functie echt willen hebben kunnen we de 
substitutieoperator /. gebruiken:

f = y[x]/.sol[[1]]
f/.x->0

Als we nu een oplossing met een bepaalde beginwaarde (hier: -2) willen hebben,
dan kunnen we natuurlijk de constante C[1] zo bepalen dat de beginwaarde gelijk wordt aan -2. Handiger is echter om Mathematica deze voorwaarde aan het
begin al mee te geven. Mathematica ziet de voorwaarde y[0] = -2 als een
extra vergelijking: we hebben nu dus twee vergelijkingen voor de onbekende
functie y[x], de differentiaalvergelijking en de beginvoorwaarde. We geven deze 
vergelijkingen net als boven aan DSolve, maar omdat we er nu twee hebben
moeten we er accolades omheen zetten.

exact = DSolve[ { y'[x] == E^(-x/3) y[x] + 1, y[0] == -2},
y[x], x ]

Wederom zien we de onbekende functie ExpIntegralEi. Bij het tekenen van de
oplossing hindert dit niet. De oplossing is weer gegeven in de vorm van een 
substitutieregel, we krijgen dus de waarde van de oplossing door de 
substitutieoperator /. te gebruiken:

plot2 = Plot[ y[x]/.exact[[1]], {x,0,20} ]

We kunnen de oplossing en het lijnelementenveld ook in één plot tekenen (de optie ImageSize kan aangepast worden om de tekening
groter, danwel kleiner te maken):

Show[ plot1, plot2, ImageSize -> 300]