Ingemar Nåsell, Institutionen för Matematik, KTH, 001102.

KURS I DIFFERENTIALEKVATIONER (DIFF OCH TRANS II, DEL 2).

OBS! Denna labb är skriven för Maple 6.

DATORLABORATION 1: Linjär diffekvation av andra ordningen med konstanta koefficienter och periodiskt högerled. Tvingade svängningar, numerisk bestämning av amplitud och fasvinkel hos periodisk lösning. Resonans.

Denna labb handlar om ett mekaniskt system som består av en massa och en fjäder med dämpning och som påverkas av en yttre periodisk kraft. Differentialekvationen för systemets avvikelse u från jämviktsläget kan skrivas

(1) .

Den oberoende variabeln står för den oskalade tiden. Ekvationen har fem parametrar: massan , dämpningskonstanten , fjäderkonstanten , yttre kraftens amplitud , och yttre kraftens vinkelfrekvens .

Med hjälp av dimensionsanalys och skalning kan man införa dimensionslösa parametrar, dimensionslös tid t och en dimensionslös tillståndsvariabel x. Efter skalningen kan diffekvationen skrivas

(2) .

De fem parametrarna har alltså reducerats till två, nämligen den dimensionslösa dämpningskonstanten och den likaledes dimensionslösa vinkelfrekvensen för den yttre kraften.

Varje lösning till denna diffekvation kan skrivas som summan av två termer, där den ena är en periodisk partikulärlösning till den inhomogena ekvationen och den andra är en transient (om > 0) lösning till den homogena ekvationen. Den periodiska lösningen är en funktion av tiden t som beror på de två parametrarna och . Amplituden och fasen hos den periodiska lösningen beror på de två parametrarna och . Den periodiska lösningen kan skrivas i formen

(3) .

Härledningen tillåter oss att bestämma amplituden och fasvinkeln som funktioner av de två parametrarna och . Uttrycket för amplituden kan skrivas

(4) ,

medan fasvinkeln satisfierar

(5) , 0 < .

Dessutom kan vi härleda ett villkor på den skalade dämpningen som garanterar att amplitudresonans kan inträffa, och vi kan visa hur resonansvinkelfrekvensen kan bestämmas som funktion av den skalade dämpningen .

Motsvarande resultat har härletts i kursboken direkt från den oskalade ekvationen (1). Vårt angreppssätt med skalning har två stora fördelar: Det ger både en ökad fysisk förståelse för problemets struktur och en förenkling av den matematiska analysen.

I denna labb studerar vi den periodiska lösningens amplitud och fas, både numeriskt och teoretiskt. Den numeriska studien består i att vi låter Maple lösa diffekvationen för givna initialvärden. Vi studerar sedan lösningen för så stora t-värden att transienten har klingat ut. Här bestämmer vi sedan amplitud och fas numeriskt och jämför dem med de teoretiskt härledda uttrycken (4) och (5).

Vi studerar differentialekvationen (2) med parametervärdena = 0.1 och = 0.9, och med initialvärdena x(0) = 1, x'(0) = 2.5. Ekvationen ges namnet ekv.

> restart;

> ekv:=diff(x(t),t$2)+2*alpha*diff(x(t),t)+x(t) = cos(beta*t);
alpha:=0.1;
beta:=0.9;

Lösningen till det ovan beskrivna initialvärdesproblemet kallas solution:

> solution:=dsolve( {ekv, x(0)=1, D(x)(0)=2.5}, x(t) );

Vi ser att solution är en ekvation (i Maples mening) med ett vänsterled och ett högerled. Vi är intresserade av uttrycket i högerledet. Det går att använda Maple för att förenkla detta uttryck, men detta behövs inte för våra numeriska undersökningar. Vi plockar ut högerledet ur uttrycket solution med kommandot rhs (right-hand side). Avsluta kommandot med ett kolon om du inte vill se uttrycket en gång till. Eftersom namnet "lösning" innehåller den svenska bokstaven ö måste namnet omslutas av bakåtlutande apostrofer.

> `lösning`:=rhs(solution):

Vi plottar lösningen över ett lämpligt tidsintervall:

> plot( `lösning`, t=0..50);

Om lösningen vore periodisk så skulle vi kunna bestämma perioden som tidsskillnaden mellan två successiva nollgenomgångar med negativ (eller positiv) derivata. Från figuren kan vi avläsa sådana t-värden med dålig noggrannhet, men Maple kan användas för att numeriskt beräkna dem med god noggrannhet. Vi använder Maple-kommandot fsolve . Detta behöver information om ett intervall där lösningen ligger. Vi väljer de två längst till höger i figuren liggande nollgenomgångarna med negativ derivata. Från figuren ser vi att dessa ligger i intervallen (36,40) och (42,46). De motsvarande t-värdena benämns t1 och t2:

> t1:=fsolve(`lösning`, t, 36..40);
t2:=fsolve(`lösning`, t, 42..46);

Skillnaden mellan dessa t-värden är perioden (om funktionen är periodisk):

> t2-t1;

Perioden hos den givna ekvationens högerled är

> period:=evalf(2*Pi/beta);

Den observerade perioden är större än perioden hos ekvationens högerled. Orsaken är att transienten inte har klingat ut än. Vi gör vår mätning för ett högre t-värde och använder periodmätningen som en koll på att t-värdet är tillräckligt högt. Man kan få experimentera med flera t-intervall innan man blir nöjd.

> plot(`lösning`, t=200..220);

Som ovan använder vi figuren för att bestämma två intervall inom vilka vi har två successiva nollgenomgångar med negativ derivata. Dessa intervall matas till fsolve nedan för att bestämma t1 och t2. Vi beräknar sedan den observerade perioden t2 - t1. I fjärde steget beräknar vi sedan kvoten mellan den numeriskt bestämda perioden och högerledets teoretiskt beräknade period.

> t1:=fsolve(`lösning`, t, 210..214);
t2:=fsolve(`lösning`, t, 218..220);
t2-t1;
kvot:=(t2-t1)/period;

Om den beräknade kvoten hade avvikit från 1 med mer än några enheter i den näst sista siffran hade detta varit en indikation på att vi borde välja vårt t-intervall längre till höger.

Vi bestämmer nu amplituden numeriskt genom att först beräkna ett t-värde där lösningen har maximum, och sedan beräkna lösningens värde vid detta t-värde. Vi börjar med att bestämma derivatan av lösningen med avseende på t.

> derivata:=diff(`lösning`,t);

Figuren ovan visar att lösningen har ett maximum i intervallet mellan 215 och 220. Vi bestämmer först t-värdet för detta maximum, d.v.s. ett t-värde i detta intervall där derivatan är lika med noll:

> tmax:=fsolve(derivata, t, 215..220);

Sedan bestämmer vi värdet på lösningen när t antar detta värde:

> xmax:=evalf( subs(t=tmax, `lösning`) );

Hur lågt är det minimum som figuren visar att lösningen har i intervallet från 212 till 216? Vi bestämmer först ett t-värde i detta intervall där lösningens derivata är lika med noll:

> tmin:=fsolve(derivata, t, 212..216);

Lösningens värde för detta t-värde är:

> xmin:=evalf( subs(t=tmin, `lösning`) );

Om absolutbeloppen av xmax och xmin hade avvikit med mera än några enheter i näst sista siffran så hade detta varit en indikation på att vi borde söka ett t-intervall med högre t-värden.

Vi jämför nu det numeriskt beräknade värdet på amplituden med det teoretiska uttrycket för amplituden, givet i (4) ovan.

> K:=1/sqrt((1-beta^2)^2+(2*alpha*beta)^2);

Överensstämmelsen mellan det numeriskt bestämda värdet på amplituden och det teoretiska värdet accepteras.

Det sista steget är att skatta fasvinkeln . Vi ritar lösningen och ekvationens högerled i samma figur:

> plot( {`lösning`, cos(beta*t)}, t=200..220);

Lösningen ges av kurvan med den större amplituden. Vi ser att lösningen ligger efter högerledet i fas. Vi skattar den nollgenomgång med negativ derivata för högerledet av kurvan som ligger närmast t2, som ju var lika med ungefär 219. Vi ser av figuren att den ligger i intervallet från 215 till 220:

> t3:=fsolve( cos(beta*t), t, 215..220);

Vi har nu bestämt t2 och t3 med och . Under antagandet att argumenten för cos-funktionerna vid de observerade nollgenomgångarna är lika har vi . Alltså kan vi bestämma fasvinkeln .

> phi:=beta*(t2-t3);

Teoretiska värdet på fasvinkeln erhålls genom att lösa , där ligger i intervallet mellan 0 och . Lösningen erhålls genom att bestämma arctan för högerledet i detta uttryck om detta högerled är > 0; eljest är lösningen + arctan för högerledet. Vi bestämmer därför först högerledet:

> `högerled`:=2*alpha*beta/(1-beta^2);

Eftersom högerledet är positivt är den teoretiskt bestämda fasvinkeln lika med arctan för detta högerled:

> arctan(`högerled`);

De åtta första siffrorna i det numeriskt bestämda värdet på fasvinkeln överensstämmer alltså med motsvarande siffror i det teoretiskt beräknade värdet på fasvinkeln ! Resultatet accepteras.