DN1240, numo08 Stefan Knutas, 881212-0056 Fredrik Båberg, 880312-0511

3B.10: Nalle-Maja gungar

Sammanfattning Detta arbete är skrivet som en del av Numeriska Metoder, Grundkurs 2. Uppgiften vi valde gick ut på att simulera när Nalle-Maja gungar på en gunga och efter att fått upp farten hoppa så långt som möjligt. Nalle-Maja gungar är en uppgift som kräver att man kan lösa differentialekvationer, i det här fallet av högre ordning, numeriskt med hjälp av MATLAB. Vi använder här Runge-Kutta samt interpolation för beräkningarna. De resultat vi kommer fram till är att vinkeln för optimalt hopp ligger runt 35°, och att man då kommer runt 2.4 meter, samt för att få ett högre resultat är det hastigheten som påverkar mest. Tillförlitligheten avgörs genom steghalvering, och vi utför även en störningsräkning för att se hur säkert programmet är för osäkerhet i indata. Våra resultat används sedan för att dra slutsatser om uppgiften, som att hastigheten är mycket viktig för att avgöra hur långt man kan hoppa.

2/6

Bakgrund Som en del i kursen DN1240, Numeriska metoder grundkurs 2, ska en större labb genomföras, samt muntligt och skriftligt redovisas. Vi har fått välja mellan ett flertal uppgifter som samtliga handlade om differentialekvationer. Detta är rapporten på den labb vi valt.

Problemställning Uppgiften var att simulera Nalle-Majas gungning och hopp, som vi delade upp i tre steg: 1. Gunga med bara dämpning 2. Ta fart i vändlägena 3. Hoppa från gungan Utöver detta skulle det längsta hoppet tas fram och hur detta kan uppnås. Givet var differentialekvationen för utslagsvinkeln och för luftfärden som då skulle omformuleras till system anpassat för MATLAB.

Följande ekvationer var givna:

d

2

u dt

2 

k m du dt



g L

sin

u

= 0 Den dämpade svängningsrörelsen (ekv. 1)

d

2

x dt

2 =−

dx dt

 

dx dt

2   

dy dt

2  Luftfärden i X-led (ekv. 2)

d

2

dt

2

y

=−

g

− ∣

dy dt

∣  

dx dt

2   

dy dt

2  Luftfärden i Y-led (ekv. 3) där k, m, g, L och κ är konstanter. L är längden på gungan (2.0 meter), m är massan för systemet (17kg), g är gravitationskonstanten (9.81), k är en dämpningsfaktor. κ är luftmotståndskoefficienten.

Som man ser i ekvationerna för luftfärden är den uppdelad i två komponenter, en i horisontellt led samt en i vertikalt. Den vertikala (ekv. 3) använder absolutbelopp på hastigheten, då den inte kommer att motverka fallet.

Tillvägagång Det första vi gjorde var att analysera problemet, för att på det sättet avgöra vilken metod vi skulle använda. Eftersom det handlade om differentialekvationer, och vi senare i uppgiften skulle införa en diskkontinuitet, så valde vi fjärde ordningens Runge-Kutta där vi kan enkelt kan påverka värdena mellan varje beräkning. Vi började med att lösa själva dämpade gungningsrörelsen (fig. 1); genom att skriva om differentialekvationen i vektorform så att det blev av första ordningen kunde den enkelt användas 3/6

tillsammans med Runge-Kutta i MATLAB 1 .

Figur 1: Dämpningen syns tydligt på de två nedersta graferna.

När vi hade skrivit koden för den dämpade rörelsen, utökade vi den med att Nalle-Maja skulle kunna ge lite extra fart i vändlägena. Vändlägena identifierades enkelt med hjälp av tecknet på hastigheten , så var gång det skiftade så kunde diskkontinuiteten implementeras. Värdet på diskkontinuitet var inte lika enkelt att välja, då det inte fanns något angivet i lydelsen. Efter att ha analyserat rörelsen ansåg vi att en diskkontinuitet på 0.7 såg relativt naturligt ut. Därefter lade vi till ytterligare ett villkor, att bara införa diskkontinuiteten om vinkeln var under 60°, dvs. om hon kommer upp över 60° så får hon inte göra fart. Just 60° valdes inte utifrån några numeriska analyser utan är bara en uppskattning till en gräns för säker gungning.

Att sedan införa själva hoppet var en större utmaning, då vi hade två differentialekvationer under luftfärden, en för acceleration i X-led, samt en i Y led. Även systemet av dessa två ekvationer löstes på samma sätt som tidigare, genom vektorform och sedan Runge-Kutta. Denna gång behövdes inga extra ingripande under lösningen så de inbyggda lösningsmetoderna så som ode45 skulle kunna fungerat lika väl, men vi valde att fortsätta med samma. Beräkningen upprepades så länge Nalle-Maja var ovan X-axeln vilket gjorde att det sista värdet hade ett negativt Y-värde. För att korrigera detta använde vi linjär interpolation mellan de två sista punktera för att då få X-värdet precis då

markytan nåddes (y=0) 2 . Begynnelsevärdena som användes var värdena från

1 Kapitel 8 “Differentialekvationer” i Numeriska algoritmer med MATLAB,

Gerd Eriksson NADA, KTH

2 Värdet får genom att ta det sista X-värdet multiplicerat med differensen mellan de två sista punkterna och subtrahera från det sista Y-värdet.

4/6

tidigare uträkning av vinkel och hastighet vid gungningen och för varje vinkel beräknades hoppet numeriskt. Därefter jämfördes alla hopp för att se vilket som var längst och vilken vinkel det utgick från (fig. 2).

När vi arbetade med den senare delen av uppgiften, att bestämma hur långt Nalle-Maja hoppar, så hamnade vi i tolkningsproblemet varifrån hoppet skulle räknas. Antingen räknar man från gungans position vid hoppet, eller från vilopunkten för gungan. När vi löste uppgiften valde vi det senare, då vi ansåg att detta var mer intressant.

Figur 2: Optimalt hopp med begränsningen 60°, i de två nedre graferna syns nu att vinkelhastigheten och vinkeln går mot en begränsning till följd av den valda diskkontinuiteten.

Tillförlitlighet Fjärde ordningens Runge-Kutta som vi använder har ett globalt fel proportionellt mot steglängden höjd till 4, vilket ger en tillräckligt bra precision vid den numeriska beräkningen, då vi har en steglängd på 0.05.

Genom att ändra noggrannheten och beräkna om kan vi få ut hur pålitlig vår beräkning är, efter en steglängdshalvering från 0.02 till 0.01 fick vi en skillnad på ~0.1° och ~0.05 meter.

För själva resultatet anser vi att noggrannheten inte behöver vara större än ±½dm, vilket vi nått med steglängden 0.05.

Felet vid linjär interpolation är obetydligt litet då interpoleringspunkterna ligger så pass nära varandra som i detta fall.

5/6

Vi utförde även störningsräkning, genom att störa vinkeln med 1% under hoppet. Detta medförde att det längsta hoppet varierade med ~0.05m vilket vi anser vara godtagbart.

Resultat Genom de beräkningar vi gjort, har vi kommit fram till att den optimala vinkeln att hoppa från är omkring 35°, eventuellt lite senare om man gungar mycket snabbt. Under dessa förhållanden kan man uppnå hopp på omkring 2.5 meter.

Slutsatser Noggrannheten i våra resultat anser vi vara mer än tillräcklig, speciellt om man väger in de naturliga omständigheterna (en björn/människa hoppandes från en gunga). Att vi hade en begränsning i vinkeln på 60° spelade inte någon roll för beräkningen, då Nalle-Maja inte kom upp i tillräckligt hög hastighet. Det är först då man testar med högre hastigheter som detta börjar spela roll. Liknande fall är det för luftmotståndet som ökar med hastigheten, man skulle kunna tänka sig uppnå en kritisk hastighet då man skulle förlora mer energi av att åka fortare. Detta inträffar dock inte för det vi anser vara naturliga hastigheter att gunga med.

6/6

Bilaga 1 Kod för enbart gungning med dämpning

close, clear, clc % Konstanter och värden på dessa L = 2.0; k = 1.20; m = 17; g = 9.81; h = 2.5; % Höjd till grenen % Funktioner fg = @(u) [u(2), -(k/m)*u(2)-(g/L)*sin(u(1))]; % Vinkelhastighet och % acceleration vid gungning % Justera grafens egenskaper subplot(2,1,1), axis([-L L 0 h]), axis equal, title('Gungan') hold on, grid on % Startvärden tslut = 50; % Tid simuleringen ska köras dt = 0.05; % Steglängd n = tslut/dt; % Antal beräkningar y = [5*pi/36, 0]; % Vinkel och vinkelhastighet rikt = -1; % Gungar åt vänster från start phi = [y(1), zeros(1, n)]; % Startvinkel samt allokering phiprim = [y(2), zeros(1, n)]; % Vinkelhastighet samt allokering for j = 1:n+1 % RK4 xL = L*sin(y(1)); yL = -L*cos(y(1))+h; plot([0 xL], [h yL], 'b-', xL, yL, 'ro') f1 = fg(y); f2 = fg(y+ dt*f1/2); f3 = fg(y+ dt*f2/2); f4 = fg(y+ dt*f3); y = y + dt*(f1 + 2*(f2+f3) + f4)/6; phi(j) = y(1); phiprim(j) = y(2); xL = L*sin(y(1)); yL = -L*cos(y(1))+h; plot([0 xL], [h yL], 'y-', xL, yL, 'go') pause(dt) % Kommentera bort för att resultatet snabbt end t = 0:dt:tslut; % Plot-variabel subplot(2,2,3), plot(t, phi, t, phiprim, ':') xlabel('Tiden'), ylabel('Vinkel, vinkelhastighet [r/s] (prickad)') subplot(2,2,4), plot(phi, phiprim) xlabel('Vinkeln \phi'), ylabel('Vinkelhastighet [r/s]') I/I

Bilaga 2 Kod för gungning och hopp

close, clear, clc % Konstanter och värden på dessa L = 2.0; k = 1.20; m = 17; g = 9.81; h = 2.5; % Höjd till grenen kappa = 0.15; % Funktioner fg = @(u) [u(2), -(k/m)*u(2)-(g/L)*sin(u(1))]; % Vinkelhastighet och acceleration vid gungning fh = @(u) [u(2), -kappa*u(2)*sqrt(u(2)^2+u(4)^2), ... % Hastighet och % acceleration vid hoppet, X-led u(4), -g-kappa*abs(u(4))*sqrt(u(2)^2+u(4)^2)]; % Y-led % Justera grafens egenskaper subplot(2,1,1), axis([-L L 0 h]), axis equal, title('Gungan') hold on, grid on % Startvärden tslut = 50; dt = 0.05; n = tslut/dt; y = [5*pi/36, 0]; % Tid simuleringen ska köras % Steglängd % Antal beräkningar % Vinkel och vinkelhastighet rikt = -1; % Gungar åt vänster från start phi = [y(1), zeros(1, n)]; % Startvinkel samt allokering phiprim = [y(2), zeros(1, n)]; % Vinkelhastighet samt allokering for j = 1:n+1 % RK4 xL = L*sin(y(1)); yL = -L*cos(y(1))+h; plot([0 xL], [h yL], 'b-', xL, yL, 'ro') f1 = fg(y); f2 = fg(y+ dt*f1/2); f3 = fg(y+ dt*f2/2); f4 = fg(y+ dt*f3); y = y + dt*(f1 + 2*(f2+f3) + f4)/6; phi(j) = y(1); phiprim(j) = y(2); xL = L*sin(y(1)); yL = -L*cos(y(1))+h; plot([0 xL], [h yL], 'y-', xL, yL, 'go') if sign(y(2)) ~= rikt % Om hastigheten har byt tecken ska det adderas % hastighet if abs(phi(j)) < pi/3 % för obegränsad gungning disp('Tar fart'); % men bara om vinkel är under 60 grader. % Denna sats (med tillhörande end) tas bort y(2) = y(2)+0.7*sign(y(2)); % Får fart på 0.7

% Används för att kunna se skillnad på när % fart inte tas vid en vändning end rikt = sign(y(2)); end I/II

pause(dt) % Kommentera bort för att resultatet snabbt end t = 0:dt:tslut; % Plot-variabel subplot(2,2,3), plot(t, phi, t, phiprim, ':') xlabel('Tiden'), ylabel('Vinkel, vinkelhastighet [r/s] (prickad)') subplot(2,2,4), plot(phi, phiprim) xlabel('Vinkeln \phi'), ylabel('Vinkelhastighet [r/s]') % Simulera längd för hopp xMax = 0; % Längsta hoppet dt = 0.05; for j = 1:length(phi) % Beräkna för varje vinkel med respektive hastighet Vy = sin(phi(j))*phiprim(j)*L; % Hastighet i Y-led Vx = cos(phi(j))*phiprim(j)*L; % Hastighet i X-led Phopp = [L*sin(phi(j)), -L*cos(phi(j))+h]; % Positionen vid hoppet w = [Phopp(1), Vx, Phopp(2), Vy]; % Position och hastighet wTemp = []; while w(3) > 0 % RK4 % För given vinkel % beräkna hoppet så länge w(3)(y) är över 0 f1 = fh(w); f2 = fh(w+ dt*f1/2); f3 = fh(w+ dt*f2/2); f4 = fh(w+ dt*f3); w = w + dt* (f1 + 2*(f2+f3) + f4)/6; wTemp = [wTemp; w]; % Spara punkterna så att de kan % ritas ut efter end % Interpolerar värdet (linjärt), för att få x då y=0 xTemp = wTemp(end,1) - wTemp(end,3)*(wTemp(end,1)-wTemp(end-1,1)) / (wTemp(end,3) wTemp(end-1,3)); if abs(xTemp) > abs(xMax); % Kontrollerar om det är längsta % hoppet hittills optimalPhi = phi(j); wMax = wTemp; xMax = xTemp; end end optimalPhiGrad = optimalPhi*180/pi; disp(strcat('Längsta hoppets nerslag [m]:', num2str(xMax))) disp(strcat('Längsta hoppets vinkel [grader]:', num2str(optimalPhiGrad))) % Rita ut längsta hoppet och skriv in värdena subplot(2,1,1), plot(wMax(:,1), wMax(:,3), 'k:') text(xMax, 0.4, strcat(num2str(xMax), 'm')) text(wMax(1,1), 0.4, strcat(num2str(optimalPhiGrad), '\circ')) II/II