1
Numerieke Analyse, Gewone differentiaalvergelijkingen. January 30, 2002
Een ski¨er wil zijn sprong van een springschans simuleren om de kritische parameters in zijn sprong te kunnen optimaliseren. Het hoogteverschil van de schans is 50 meter en de vorm van de schans is een vloeiende kromme hoogte = ϕ(x) met een beginhelling van 45 graden, met ϕ00 ≥ 0 en met een horizontaal einde. Het landingsgebied heeft een constante helling van 30 graden en het hoogteverschil tussen het eind van de schans en het begin van het landingsgebied is 2 meter. De lengte van een sprong wordt gemeten langs een horizontale lijn die begint loodrecht onder het eind van de schans.
Leid differentiaalvergelijkingen af voor de positie in IR2 met x de horizontale en y de verticale component, als je voor de wrijvingskracht een model neemt met fw = a v + b v2 waarin v de snelheid is. De lineaire term modelleert de baanweerstand en de kwadratische de luchtweerstand. Kies de parameters a en b zo, dat de energiedissipatie (bij het glijden op de schans) bij een snelheid van 36 km/u (10 m/s) 150 watt is en dat de bijdragen aan de energiedissipatie van de lineaire en de kwadratische term (ongeveer) gelijk zijn bij een snelheid van 3 m/s. Tijdens de vrije val is a gelijk aan nul. Tijdens de vrije val kunnen we nog een extra wrijving modelleren ten gevolge van de luchtstroom tegen de onderzijde van de ski’s, die een hoek α(t) met de voortbewegingsrichting maken. Deze kracht staat loodrecht op het vlak van de ski’s en is evenredig met het oppervlak van de ski’s maal de sinus van de hoek met de voortbewegingsrichting. Het zou te ver voeren om dit als optimaal besturingsprobleem te behandelen. Kies er daarom voor dat α constant is t.o.v. de horizon. Gebruik van de wet van Newton levert twee tweedeorde vergelijkingen voor de x- en de y-coordinaat van de positie van de skier; zet dit om in een stelsel van vier vergelijkingen van eerste orde voor u = (x , ˙x , y , ˙y).
Als de wrijving nul is, is de snelheid bij het verlaten van de schans eenvoudig te bepalen en is ook de baan daarna en het punt van neerkomen uit te rekenen. De uitkomst is dan ook niet afhankelijk van de massa van de ski¨er (we beschouwen deze als een “massapunt” in het gehele model). Met het geval dat de wrijving nul is kunnen we testen of onze methode een goede uitkomst geeft en of de afbreekfout op de goede wijze afneemt bij het halveren van de stapgrootte.
Als de wrijving niet nul is, is er geen exacte oplossing en zullen we de lengte van de sprong met een numerieke methode moeten benaderen. We moeten eerst rekenen met het schansmodel tot het tijdstip dat de ski¨er de schans verlaat. Vervolgens moeten we rekenen tot het moment dat de ski¨er op de grond komt. In beide gevallen moeten we integreren tot voorbij het moment van loslaten resp. neerkomen.
Het is de bedoeling een nauwkeurige berekening te maken van de baan en in het bijzonder van het punt van neerkomen.
Groep 1 gebruikt de methode van Heun met constante stapgrootte h. Geef ook een bewijs dat de lokale afbreekfout van orde h2 is.
Groep 2 gebruikt Adams-Bashforth orde 3 predictor en Adams Moulton orde 4 corrector met constante stapgrootte. Geef ook een bewijs dat de lokale afbreekfout van orde h4 is.
We ontwikkelen en checken de methode in een aantal stappen, te beginnen met het glijden tot het tijdstip van loslaten t∗. We kiezen (0,0) in het xy-vlak als de positie van het eind van de schans en we laten de beweging van links naar rechts verlopen; d.w.z. de schans heeft x < 0 en y > 0 en de helling waarop de ski¨er neerkomt heeft x > 0 en y < 0. We laten de ski¨er op t = 0 bovenaan de schans beginnen.
1. Het testen van de integratiemethode.
We kiezen een initi¨ele stapgrootte h0 en we berekenen het laatste tijdstip bt = k hb 0 waarop de benadering xk van x(bt) (berekend met deze stapgrootte) nog positief is. Vervolgens bereken je nauwkeuriger benaderingen van x(bt) voor stapgrootten h = h0/2 , h = h0/4 , · · ·.
Groep 1 doet hierop h2-extrapolatie en (indien nodig) h3-extrapolatie en checkt of de fout in de zo verkregen rij extrapolaties inderdaad van orde O(h4) is.
Groep 2 hoeft alleen maar na te gaan of de (geschatte) fout in de verkregen rij benaderingen van orde O(h4) is, zoals de theorie zegt.
2. Het benaderen van t∗.
Als we het stelsel integreren met stapgrootte h vinden we een k zodat voor de benadering van x op tijstip tk:= kh geldt xk < 0 en xk+1 > 0 . Het moment van loslaten heeft dus plaatsgevonden voor een th∗ ∈ (kh , kh + h).
Oefeningen Numerieke Analyse, ODE 2
Groep 1 berekent t∗ door inverse lineaire interpolatie; d.w.z. neem xk en xk+1 als abscissae en kh en kh + h als bijbehorende functiewaarden en definieer th∗ als de waarde van het lineaire interpolatiepolynoom in x = 0. Laat zien dat de afbreekfout in th∗ van orde O(h2) is evenals de afbreekfout in xk.
Groep 2 berekent ook het virtuele punt xk+2, alsof de schans horizontaal blijft doorlopen voor positieve x; aangezien de ski¨er een positieve horizontale snelheid heeft, zal gelden xk−1< xk<
xk+1< xk+2 en kunnen we th∗ berekenen door inverse kubische interpolatie(dus de waarde van het derdegraads interpolatiepolynoom in x = 0 door de punten (xj, jh)). Laat zien, dat dit een afbreekfout in th∗ van orde O(h4) garandeert.
Waarom is interpolatie op xk−1< xk< xk+1< xk+2 beter dan op xk−2< xk−1< xk< xk+1? Doe deze berekening voor de stapgrootten h = h0, h = h0/2 , h = h0/4 ,· · · en onderzoek of h2- extrapolatie resp. h4-extrapolatie op de gevonden rij benaderingen mogelijk is; d.w.z. onderzoek of voor de afbreekfout geldt th∗ = t∗ + chp + O(hp+1) (met p=2 resp. p=4) of niet meer dan
| th∗ − t∗| ≤ chp en probeer vervolgens dit resultaat aannemelijk te maken.
3. De vrije val en het neerkomen.
Integreer vanaf th∗ met dezelfde stapgrootte verder en bepaal op dezelfde manier als in het vorige punt het tijdstip theindvan neerkomen en de bijbehorende positie (i.h.b. dus de lengte van de sprong).
Verifieer of de afbreekfout nog steeds van orde O(hp) is.
4. Maak voor de beste waarde van h een grafiek van de baan in het xy− vlak en een grafiek met de componenten van u als functie van de tijd (eventueel via vier subplots, zie “help subplot” in Matlab).
• Ter vergelijking kun je desgewenst de standaard-matlabroutines ode23 en ode45 gebruiken.
• Een voorbeeld van een extrapolatietabel vind je in de syllabus in de paragraaf over Romberg- integratie.
• Schrijf voor de overzichtelijkheid aparte functies voor de integratieroutine en voor de (vector-)functie van de differentiaalvergelijking.
• Om een functie (gedefinieerd door een m-file) in een andere functie te kunnen oproepen heb je het mechanisme uit het bijgaande voorbeeld nodig. Je importeert de naam van de functie (als een variabele van het type string) en je maakt een string door concatenatie van de functienaam met een string van actuele parameters. Matlab kan deze string dan evalueren alsof het een statement is via de functie “eval”.
function y=integraal(fun,a,b,n)
% oproep I=integraal(’mijnfun’,0,1,23)
% voor het integreren van de functie met naam "mijnfun".
% Deze functie kan een standaardfunctie zijn in matlab
% zoals "sin" of "atan" of "exp" of een zelfgeschreven
% m-file met naam "mijnfun.m", zie beneden
% de functie moet zo geschreven zijn dat deze een vector als
% parameter kan meekrijgen en een even grote vector als resultaat
% aflevert h=(b-a)/n;
x=a+(0:n)*h;
funexpr=[fun,’(x)’]; % concatenatie van funtienaam met de actuele parameters z=eval(funexpr); % evaluatie van de string als matlab-opdracht
y=(z(1)+z(n+1)+2*sum(z(2:n)))*h/2;
% Voorbeeld van een functiefile:
%% function w=mijnfun(x);
%% w=1./(x.^2+1)
% Als je m-file-functie alleen kan rekenen met een scalaire
% parameter, dan moet je de laatste drie regels van de integraal-m-file
% vervangen door:
%% z=zeros(size(x));
%% for k=1:n+1,funexpr=[fun,’(x(k))’];z(k)=eval(funexpr);end
%% y=(z(1)+z(n+1)+2*sum(z(2:n)))*h/2;
