Gewone differentiaalvergelijkingen: Lijnelementveld en oplossingskrommen met MATLAB
Tekenen van een lijnelementenveld met de quiver functie in MATLAB
Stel dat we een differentiaalvergelijking van de vorm hebben voor een zekere functie en dat we het lijnelementenveld willen tekenen m.b.v. MATLAB. Op elk punt kunnen we dan de helling van een oplossing berekenen en in het -vlak een kort lijnstuk met deze helling tekenen. Om dit op een regelmatig rooster te doen zijn de functies meshgrid en quiver handig.
Basissyntax voor lijnelementenveld incl. voorbeeld De basissyntax voor het tekenen van een lijnelementenveld in MATLAB is als volgt:
[t,y] = meshgrid(tmin:tstep:tmax, ymin:ystep:ymax);
dy = phi(t,y);
dt = ones(size(dy));
figure
quiver(t, y, dt, dy)Bij de coderegel waarin dt gelijk gesteld wordt aan allemaal eentjes, kun je denken aan het schrijven van de differentiaalvergelijking als Een concreet voorbeeld met ziet er als volgt uit:
>> [t,y] = meshgrid(0:0.12:2, 0:0.1:2);
>> dy = 1 - 2*t.*y;
>> dt = ones(size(dy));
>> figure
>> quiver(t, y, dt, dy);
Het resultaat is een beetje teleurstellend omdat de functie quiver automatisch de pijltjes schaalt zodat ze altijd in beeld zijn, maar dat de pijlen wel sterk afwijkende lengtes hebben. Ook is er veel onnodige witruimte boven en aan de rechterkant.
Merk op dat we eigenlijk een vectorveld hebben getekend (de verstekwaarde van quiver) en niet een lijnelementenveld. Ook is het vertrekpunt van een pijl in het roosterpunt en leggen we niet het midden van een vector (of lijnelement) op het bijpassende roosterpunt. Via de optie ShowArrowHead (verstekwaarde 'on') kunnen we de vectoren vervangen door lijnstukjes, maar we zullen in het vervolg hiet niet om malen omdat de pijltjes ook de richting van een oplossingskromme aangeven.
De grootste verbetering in het lijnelementveld krijgen we door zelf de pijlen te schalen zodat ze allemaal gelijke lengte hebben.
Lijnelementenveld met schaling van lijnelementen Schaling van de vectoren naar lengte 1 kan in het sjabloon via de commando's
>> L = sqrt(dt.^2 + dy.^2)
>> dyu = dy./L;
>> dtu = dt./L;
We weten nu zeker dat de vectoren in allemaal lengte 1 hebben. Met de optie axis tight halen we overtollige witruimte in het lijnelementenveld weg. Ons eerdere voorbeeld wordt nu de volgende verbeterde figuur:
>> clear;
>> [t,y] = meshgrid(0:0.1:2, 0:0.1:2);
>> dy = 1 - 2*t.*y;
>> dt = ones(size(dy));
>> L = sqrt(dt.^2 + dy.^2);
>> dyu = dy./L;
>> dtu = dt./L;
>> figure
>> quiver(t, y, dtu, dyu), axis tight
Lijnelementenveld met verschillende opties Je kunt via opties de tekening van het lijnelementenvled naar wens aanpassen. We geven een voorbeeld waarin we de lengte van de pijltjes verkleint hebben met een schaalfactor gelijk aan 0.5, de kleur rood hebben gekozen en label plus een titel hebben toegevoegd.
>> clear;
>> [t,y] = meshgrid(0:0.1:2, 0:0.1:2);
>> dy = 1 - 2*t.*y;
>> dt = ones(size(dy));
>> L = sqrt(dt.^2 + dy.^2);
>> dyu = dy./L;
>> dtu = dt./L;
>> figure
>> quiver(t, y, dtu, dyu, 0.5, 'r'), axis tight
>> xlabel 't', ylabel 'y';
>> title("Lijnelementenveld bij y' = 1 - 2ty",'Color','Red')
