Kontinuerliga system, Datorövning 2

(1)

V˚arterminen 2003

Kontinuerliga system, Dator¨ ovning 2

Inledning

Ett modernt datoralgebrasystem har som huvudfunktion att göra symboliska beräkning- ar, i motsats till numeriska. Det kan utföra algebraiska manipulationer och förenklingar, lösa ekvationssystem, integrera och derivera symboliskt och lösa differentialekvationer.

Dessutom kan de flesta s˚adana system ocks˚a utföra numeriska beräkningar och har kraft- fulla grafiska funktioner. Behöver man göra mera omfattande numeriska beräkningar använder man sig dock hellre av Matlab eller specialskrivna numeriska program.

Maple är ett av de ledande datoralgebrasystemen. N˚agra av de starkaste konkurren- terna är Mathematica, Macsyma och Reduce. LTH har en generell licens för Maple, och systemet bör därför finnas tillgängligt p˚a alla institutioner och p˚a skolans elevdatorer.

Det kan k¨oras p˚a de flesta vanliga operativsystem.

Syftet med denna datorlaboration ¨ar att Du skall bli bekant med Maple och se hur det kan anv¨andas

• f¨or symboliska och numeriska ber¨akningar,

• f¨or visualisering,

• i specialanv¨andningar, genom ett antal utbyggnadspaket, s˚a kallade ’packages’.

Hinner du inte med alla uppgifterna under den handledda ¨ovningen, s˚a genomf¨or de

˚aterst˚aende momenten p˚a egen hand, antingen p˚a LTH:s datorer eller p˚a egen.

Vi kommer i denna övning att arbeta med en del exempel fr˚an övningar och föreläs- ningar. Tag därför med läroboken och övningshäftet till laborationen. Läs noga igenom den allmänna Maple-introduktionen i avsnitt 1, s˚a noga att du bara behöver ägna högst en halvtimme av laborationstiden ˚at detta. I avsnitt 2 och 3 finns n˚agra förberedelseupp- gifter. Gör dessa innan laborationen. Skriv in dina svar p˚a sista sidan och var beredd att motivera svaren under laborationen.

1 Om Maple

Allm¨ant

Maple är ett mycket m˚angsidigt program, och man kan ägna mycket tid ˚at att utforska dess användningsmöjligheter. Dess m˚angsidighet gör det samtidigt sv˚arare att använda

¨an tex Matlab.

I denna övning kommer ganska mycket utrymme att ägnas ˚at allmänna förberedel- ser, innan vi kommer till de speciella användningarna inom Kontinuerliga system. Där kommer vi att använda Maple som en oerhört kraftfull räknedosa.

Liksom de flesta programpaket uppgraderas Maple relativt ofta. Versionerna skiljer sig n˚agot ˚at, och anvisningarna nedan ¨ar gjorda f¨or Maple6. Mer information om Maple

än vad som ryms i denna korta introduktion kan man finna i det inbyggda hjälpsystemet samt i en stort antal skrifter och böcker. P˚a institutionen har utgivits en utförligare introduktion (mest avsedd för envariabelanalys) Andersson–Nilsson, Maplehandboken (1995). Den beskriver version MapleVR3 och finns att hämta p˚a

(2)

http://www.maths.lth.se/matematiklu/personal/goran_a/Maple/mh.html.

En modernare version har f¨orfattats av C.G. Werner: Maplehandboken i Lund 2002.

Häftet kan köpas p˚a studerandeexpeditionen vid Matematikavdelningen (NF) rum MH 456. En handledning för MapleVR5 är Heal, Hansen & Rickard, Maple V Learning Guide, Springer 1998.

Starta Maple p˚a UNIX-system genom att ge kommandot maple -x, eller, p˚a vissa maskiner, xmaple &. P˚a Windows-maskiner brukar finnas en Maple-ikon. Placera Maplefönstret, som strax kommer upp, p˚a lämpligt ställe. Man g˚ar ur Maple genom att välja Exit i File-menyn. Behöver man avbryta n˚agon process kan man använda en stop-knapp p˚a menyraden.

Hj¨alp

Maple har ett omfattande inbyggt hjälpsystem. Om man klickar p˚a Help i menyraden, s˚a f˚ar man tillg˚ang till hjälpfunktioner. Dessa fungerar p˚a lite olika sätt, beroende p˚a vilken Mapleversion man kör, men brukar vara lätta att använda.

Plocka fram hjälpen för exp och constant. Observera i det senare fallet särskilt de symboliska konstanterna I, Pi och infinity. Dessa är reserverade namn, som inte f˚ar användas till annat. Observera att Maple skiljer p˚a stora och sm˚a bokstäver. Talet e skrivs exp(1).

Man kan ocks˚a f˚a hjälp genom att skriva ett sökord föreg˚anget av ett fr˚agetecken, t ex ?exp, ?constant. I slutet p˚a varje hjälptext finns exempel och hänvisningar till andra kommandon som kan vara av nytta.

Grundl¨aggande kommandon och aritmetiska ber¨akningar

I Maplef¨onstret hittar du en ’prompt’ >, som betyder att Maple v¨antar p˚a ett kommando.

Skriv in 7+5; och tryck p˚a Enter, s˚a utf¨ors summationen.

Observera att varje kommando till Maple m˚aste avslutas med ett semikolon ;, utan detta händer ingenting. Detta gör att man kan sl˚a in l˚anga formler som inte f˚ar plats p˚a en rad. Skriver man kolon : s˚a utförs operationen, men resultatet kommer inte upp p˚a skärmen.

Kontrollera att Maple kan fungera som en vanlig r¨aknedosa. Prova i tur och ordning 2-3*5; 1/2+1/3; sqrt(4)*sqrt(3);

Observera att Maple ger exakta svar, inte n¨armev¨arden.

Närmevärden erh˚alles genom kommandot evalf;, som st˚ar för ’evaluate using flo- ating point arithmetic’. För att f˚a ett närmevärde för den senast utförda beräkningen (2√

3 i v˚art fall) kan man skriva

evalf(%);

Man kan naturligtvis ocks˚a skriva evalf(sqrt(4)*sqrt(3));. P˚a samma sätt kan resultatet av den näst senaste beräkningen ˚aberopas genom %%, etc.

En annan praktisk detalj i Maple är möjligheten att g˚a tillbaka till gamla kommandon p˚a skärmen med pil-upp och pil-ner tangenterna, ändra i kommandona och utföra dem p˚a nytt. Man kan ocks˚a använda musen för att snabbt flytta sej till olika delar av skärmen.

(3)

Testa genom att g˚a tillbaka till kommandot 1/2+1/3; och gör en annan br˚akberäkning genom att ändra siffror.

När det gäller kommandot evalf s˚a kan man ur hjälpfunktionen, som n˚as med

?evalf, utläsa att man genom ett andra argument kan ange med hur m˚anga siffror man vill ha svaret. Använd detta för att f˚a ett värde p˚a π med 1000 decimaler, evalf(Pi,1000);.

Variabler

Maple kan inte bara räkna med tal utan ocks˚a med variabler och med funktioner. Detta gör att ett Maplesystem blir större och ofta mer invecklat att programmera än ett vanligt programmeringspr˚ak, men ocks˚a oerhört mycket mera kraftfullt.

Ge kommandot (x+1)^3;. Maple svarar med samma sak. Som svar p˚a expand(%);

s˚a utvecklar Maple uttrycket enligt binomialteoremet. Prova ¨aven med

expand((a-b)*(a+b)); expand(x*(x+1)*(x+2)*(x+3)); expand(cos(x+y));

Som bekant vill man ofta g˚a ˚at andra h˚allet, och faktoruppdela ett givet uttryck.

Detta g¨ors med kommandot factor. Prova detta p˚a de polynomuttryck som du fick ovan.

Ett annat kommando för förenkling är simplify. Prova detta genom att skriva simplify(1/(x-1)+1/(x+1));

Maple gör inte alltid det man önskar. Ett annat kommando man kan ha nytta av för att göra omskrivningar är normal.

Tilldelningssatser

Man kan tilldela en variabel ett v¨arde, numeriskt eller symboliskt. Tilldelningssymbolen

¨ar liksom i Pascal och Simula :=. Prova med x:=2;. Kommandot x; ger nu variabelv¨ardet 2. Prova ocks˚a med (x+1)^2;.

En variabel som f˚att ett värde beh˚aller detta tills man g˚ar ur Maple eller ger den ett annat värde. För att ta bort värdet kan man använda kommandot x:=’x’;, som gör att Maple tolkar bokstaven x som en variabel, betecknad x, inte som variabelns värde.

Det är lätt att glömma bort att man gett en variabel ett värde tidigare, vilket kan leda till obegripliga resultat av räkningar. Vill man ta bort alla värden p˚a variabler och helt och h˚allet börja om fr˚an början skriver man restart;.

Variabler kan förutom numeriska värden även ha Mapleuttryck som värden. Genom tilldelningskommandot

f:=(x+1)^3;

s˚a sätter vi f lika med (x + 1)³. (Blir svaret 27 s˚a tag bort det tidigare värdet fr˚an x och försök igen.) Man kan sedan räkna vidare med f och skriva t ex f^2; och expand(f^2);.

Om man vill beräkna värdet av uttrycket f^2 för t ex x=2 utan att varaktigt tilldela xett värde s˚a kan man skriva

subs(x=2,f^2);

(4)

Funktioner

Maple kan hantera inte bara analytiska uttryck, som vi sett exempel p˚a ovan, utan

även funktioner. Dessa kan definieras p˚a flera olika sätt. Det för v˚ara ändam˚al enklaste p˚aminner om beteckningen

x 7→ g(x) f¨or en funktion g. Ge kommandot

g:= x->(x+1)^3;

Nu kan funktionsvärden beräknas p˚a det sätt man är van vid. Prova t ex med g(0);, g(-1);, g(a); och g(y+z);.

Anmärkning Man m˚aste noga h˚alla isär begreppen funktion och funktionsuttryck, där det senare behandlades under rubriken ’Tilldelningssatser’ ovan. För säkerhets skull upprepar vi skillnaden:

Om vi ger ett värde till F genom tilldelningen F:=exp(x)-sin(x); s˚a är F ett uttryck med variabeln x inbyggt. Om vi i stället ger F ett värde genom F:=x->exp(x)-sin(x);

s˚a är F en funktion, där x bara används för att definiera en regel, och där man skulle kunna använda vilken annan symbol som helst. Funktioner är mycket mer flexibla, men ibland n˚agot mer sv˚arhanterliga än funktionsuttryck. Det finns lyckligtvis ett enkelt sätt att göra om ett uttryck till en funktion, nämligen genom att använda kommandot unapply. Till exempel ger kommandot F:=unapply(exp(x)-sin(x),x); samma sak som F:=x->exp(x)-sin(x);.

Maple klarar ocks˚a funktioner av flera variabler. Funktionen h(x, y) = x²− y³+ x²y³ definieras i Maple p˚a f¨oljande s¨att:

h := (x,y) -> x^2-y^3+x^2*y^3;

Vad ger h(0,0); och h(3,2);?

Derivation

Derivator ber¨aknas med kommandot diff. Anv¨andningen framg˚ar ur exemplen diff(x^3+2*x^2,x);och diff(tan(x),x);

Andraderivator ber¨aknas enligt

diff(tan(x),x,x);

P˚a samma sätt beräknas derivator av högre ordning, där antalet x anger ordningen. För detta finns ett kortare skrivsätt. Beräkna fjärdederivatan av tan(x) genom att skriva

diff(tan(x),x$4);.

Prova även att förenkla svaret med hjälp av factor(%);.

Vi kan ocks˚a derivera det tidigare definierade funktionsuttrycket f och funktionen g genom att skriva

diff(f,x); resp diff(g(x),x);

(5)

T¨ank ut vad som kommer ut ur kommandona

diff(f,y); , diff(g,x); , diff(g(y),y);

och prova om det st¨ammer.

Anm¨arkning. Maple har ocks˚a en derivationsoperator D, se ?D. Pr¨ova tex D(g);

och D(sin);.

Man kan även lätt beräkna partiella derivator. Till exempel f˚ar man

∂h

∂x och ∂h

∂y f¨or funktionen h ovan genom kommandona

diff(h(x,y),x); resp diff(h(x,y),y);

Högre derivator f˚as genom uppräkning av variabelnamnen. Pröva t ex diff(h(x,y),x,y);

Integration

Med hjälp av kommandot int(f,x); kan man beräkna en primitiv funktion till funktionen f (x) med avseende p˚a x. Integrationskonstanten f˚ar man sedan själv lägga till om det skulle behövas. Pröva

int(f,x); , int(g(x),x);

Repetera nu n˚agra enkla primitiva funktioner fr˚an envariabelanalysen: cos x (med kommandot int(cos(x),x);), e^t (skrivet exp(t)), och 1

√1 − z² med avseende p˚a re- spektive variabel. Kvadratrotfunktionen skrivs i Maple som sqrt.

Om Maple inte kan n˚agon primitiv funktion, s˚a svarar Maple genom att ge tillbaka det man stoppat in. F¨ors¨ok med int(x^x,x); s˚a f˚ar du tillbaka

Z x^xdx Anm¨arkningProva ocks˚a

Z r 1 − 2x²

1 − x² dx och Z

e⁻^x²dx .

Dessa exempel visar att Maple känner till fler primitiva funktioner än de vanliga ele- mentära funktionerna. Den första av integralerna är exempel p˚a en s˚a kallad elliptisk integral. Vissa elliptiska integraler finns i Maple, se t ex ?EllipticE. I svaret p˚a den andra integralen ing˚ar den s˚a kallade ’error function’, erf, se läroboken kapitel 4. Denna uppträder bl a i samband med värmeledning och inom statistiken.

För att beräkna bestämda integraler, till exempel Z ³

−2

(x + 1)³dx

(6)

m˚aste även gränserna anges. Detta görs i Maple genom x=-2..3, eller allmänt x=a..b för a ≤ x ≤ b. Detta har vissa likheter med Matlabs a:b, men fungerar b˚ade för diskreta och kontinuerliga variabler. Integralen ovan beräknas genom

int((x+1)^3,x=-2..3);.

Gränserna f˚ar gärna inneh˚alla variabler, försök t ex med int(f,x=-y..(z+1));.

Aven generaliserade integraler kan ber¨aknas. Prova t ex med¨

int(1/(1+x^2),x=0..infinity); och evalf(int(1/(1+x^2),x=0..infinity));

Se ¨aven vad som h¨ander med Int(1/(1+x^2),x=0..infinity);

Summation

Maple kan ocks˚a användas för att beräkna summor. Syntaxen för summation är naturlig, P_n

k=ma_k f˚as genom Maplekommandot sum(a(k),k=m..n).

Ber¨akna potenssummorna

n

X

k=1

k och

n

X

k=1

k²

med hj¨alp av kommandona sum(k,k=1..n); och sum(k^2,k=1..n);. F¨orenkla med simplify.

Aven vissa seriesummor kan Maple ber¨akna. Vad blir¨

∞

X

k=1

1 k⁴ ? Ekvationsl¨osning

Maple kan användas för att lösa ekvationer, se ?solve och ?fsolve. Prova de tre vari- anterna

solve(x^3+x^2+x-3=0,x);

fsolve(x^3+x^2+x-3=0,x);

fsolve(x^3+x^2+x-3=0,x,complex);

och läs om dem i hjälpfunktionen. Pröva ocks˚a

solve(sin(x)=1/2,x);

Observera att i det senare fallet ger Maple i grundinställning endast en lösning. Men skriv _EnvAllSolutions:=true; Lös ekvationen sin x = ¹₂ igen. Saknas n˚agra lösningar nu? För polynomekvationer ger fsolve med optionen complex alla lösningar.

Differentialekvationer

Maple kan ocks˚a lösa differentialekvationer, se ?dsolve. Pröva med att lösa y⁰⁰+ λy = 0

med kommandot

dsolve(diff(y(x),x$2)+lambda*y(x)=0, y(x));

(7)

Grafik

Maple har omfattande grafiska m¨ojligheter.

Om du i fortsättningen vill att Maples figurer hamnar i egna fönster (vilket kan vara lättare att hantera) s˚a g˚ar du in under menyn Options, väljer Plot Display och sen Window.

Med kommandot plot ritar man tv˚adimensionella figurer i Maple, t ex funktions- kurvor. Se hj¨alpen till kommandona plot. Prova med

plot(sin(x),x=-Pi..Pi);

Axlarna skalas automatiskt s˚a att bilden fyller grafikfönstret maximalt. Om man vill ha samma skala p˚a b˚ada axlarna kan man g˚a in p˚a menyraden vid Projection och välja Constrained. För att rita flera kurvor i en figur använder man [ ]-parenteser ( även {}-parenteser fungerar). Pröva

plot([sin(x),x-x^3/6],x=-Pi..Pi);

Vad ¨ar det f¨or speciellt med polynomet?

Man kan ocks˚a rita kurvor i parameterform, pr¨ova t ex

plot([cos(x),sin(x),x=-Pi..Pi],scaling=constrained);

Hur ser kurvan ut om scaling=constrained tas bort?

Antalet grafikfönster växer snabbt. Man stänger dem man inte behöver genom att g˚a in rutan i vänstra övre hörnet i grafikfönstret märkt med − och dra till Close.

Maples tredimensionella grafik har en mängd varianter som vi inte g˚ar in p˚a här. Läs om dem i ?plot3d. Funktionsytan z = xe⁻^x²⁻^y² kan ritas med

plot3d(x*exp(-x^2-y^2),x=-2..2,y=-2..2);

Man kan vrida och vända p˚a bilden genom att placera musmarkören i grafikfönstret och h˚alla vänster musknapp nere och flytta p˚a musen l˚angsamt.

Prova med hjälp av menyer de olika möjligheterna för ’axes’, ’color’ och ’style’.

Det g˚ar ocks˚a bra att beskriva ytor i parameterform. Cylindern y² + z² = 1 l¨angs x-axeln kan ritas med

plot3d([s,cos(t),sin(t)],s=-2..2,t=0..2*Pi);

V¨alj h¨ar Projection Constrained.

Maple har m˚anga fler grafikrutiner. För att f˚a tillg˚ang till dessa görs kommandona with(plots); och with(plottools); Här är plots och plottools namnet p˚a ett s˚a kallat ’package’. Det finns ett antal s˚adana i Maple, för olika användningar.

Spara och skriva ut

Man kan spara all text i Maple-fönstret genom att klicka p˚a File och Save as samt ange ett filnamn. Texten sparas d˚a i en fil med namnet filnamn.mws. När du senare vill ta fram texten igen klickar du p˚a File och Open och anger hela filnamnet, även .mws. Det g˚ar ocks˚a att samla alla kommandon i en scriptfil med ett enkelt namn, tex prov. För att köra scriptet skriver du i Maple-fönstret read prov;.

Om du vill ta ut en bild p˚a skrivare, kan du spara den i Post-Script format genom att klicka p˚a File och Print samt ange filnamn.ps.

(8)

2 Skal¨ arprodukter, normer och projektioner

Starta med att skriva restart; för att ta bort gamla värden p˚a f, g mm. Vi kommer att beräkna m˚anga skalärprodukter och det är därför lämpligt att definiera skalärprodukten som en funktion av tv˚a variabler f och g. L˚at oss kalla denna skal. Skriv

skal:=(f,g)->int(f*g*w,x=x0..x1);

Innan denna anv¨ands m˚aste man specificera viktsfunktionen w och integrationsgr¨anserna x0 och x1.

1. Förberedelser: Skriv upp formeln för skalärprodukten (f | g) i L²([0, L]). Skriv koefficienterna c_k i en ortogonalutveckling, P

kc_kϕ_k, med hj¨alp av skal¨arprodukter.

Gör en cosinusutvecklingen av en funktion f (x) p˚a intervallet [0, L], välj {ϕk(x)}^∞⁰ = {cos kπx/L}^∞⁰ och bestäm koefficientern med kommandot skal. Starta med att sätta

w:=1; x0:=0; x1:=L;

Tala d¨arefter om att k och j ¨ar heltal genom att skriva assume(k,integer, j, integer);.

Kontrollera ortogonaliteten genom att skriva

skal(cos(k*Pi*x/L),cos(j*Pi*x/L));

Pr¨ova ocks˚a fallet med j = k och speciellt j = k = 0. Definiera koefficienterna c_k i cosinusutvecklingen som

c:=k->skal(cos(k*Pi*x/L),f(x))/skal(cos(k*Pi*x/L),cos(k*x*Pi/L));

Skriv c(k); för att f˚a en integralformel för c_k. Undersök speciellt c(0); (När det finns villkor p˚a k markerar Maple det med ett ∼ efter variabeln, se ?assume.)

Fr˚aga: Hur stämmer detta med formlerna för koefficienterna i cosinusutvecklingen p˚a sidan 62 i boken? Händer n˚agot speciellt för k = 0? I s˚a fall vad och varför?

Testa svaret i ¨ovning 3.2 a) genom att v¨alja f (x) = x och L = 1. Avsluta med att

˚aterst¨alla k, j med k:=’k’; j:=’j’;.

2. Vi skall nu lösa uppgift H.22 i övningshäftet med hjälp av Maple.

Förberedelser:Skriv upp formeln för den skalärprodukt som skall användas i H.22 och ange w, x0 och x1. Skriv med hjälp av projektionsformeln upp en formel för det sökta polynomet p. Uttryck koefficienterna med skalärprodukter och tänk efter vilka integraler som skall beräknas för att lösa uppgiften.

B¨orja med att skriva in viktsfunktion w och integrationsgr¨anserna x0 och x1. Defi- niera sen

L0:=1; L1:=1-x; L2:=x^2-4*x+2; f:=exp(-alpha*x);

Plotta de tre polynomen i samma figur (se sidan 7 i denna handledning). Kontrollera med hj¨alp av skal att de ¨ar ortogonala i L2((x0, x1), w).

Prova med att beräkna tex skal(L0,f);. Maple vill inte beräkna integralen (var- för?). Man m˚aste ställa krav p˚a α för att integralen skall vara konvergent. Skriv därför

assume(alpha>-1/2);

(9)

Kan du ber¨akna skal(L0,f); nu? Ber¨akna koefficienterna c0, c1, c2 , t ex c0:=skal(L0,f)/skal(L0,L0);

Ber¨akna sedan p = c0L0+ c1L1+ c2L2 genom att skriva p:=c0*L0+c1*L1+c2*L2;

Plotta funktionen och approximationen f¨or n˚agra olika v¨arden p˚a α, tex

plot(subs(alpha=1,[f,p]),x=0..10, color=[red,blue], thickness=2);

Fr˚aga: I vilken mening ¨ar p en bra approximation av e⁻^αx? F¨or att undvika problem senare, ta bort villkoret p˚a α med

alpha:=’alpha’;

3 Ortogonalpolynom

I Maple finns flera typer av ortogonalpolynom i ett s¨arskilt package. F¨or att f˚a tillg˚ang till dem, skriv

with(orthopoly);

3. Kontrollera med hjälpfunktionen att de ortogonalpolynom som finns nämnda i läroboken sid 271–273 ocks˚a finns i Maple. För att f˚a ett explicit uttryck för t ex Legendrepolynomet P4 skriver man P(4,x); och för att f˚a Laguerrepolynymet L2 skriver man L(2,x);.

Jämför Maples Laguerrepolynom med de som finns i övn H.22? Finns n˚agon skillnad?

Ber¨akna n˚agra skal¨arprodukter med skal, t ex skal(L(2,x),L(2,x));. Kommentar?

Plotta n˚agra ortogonalpolynom, t ex

plot([P(3,x),P(4,x),P(5,x)],x=-1..1);

Detta kan ocks˚a skrivas plot([P(n,x)$n=3..5],x=-1..1); Kan man med hjälp av fi- gurerna avgöra om det är rimligt att polynomen är ortogonala? Detta är lättare att bedöma om man tittar p˚a produkterna P4P5 och P3P5, i intervallet [−1, 1]. Undersök ortogonaliteten med hjälp av skal. (Tänk p˚a att först ange rätt viktsfunktion och integ- rationsgränser.)

4. Utveckla funktionen f = sin(10x) med avseende p˚a Legendrepolynom ¨over intervallet

−1 ≤ x ≤ 1. S¨att

f:=sin(10*x);

H¨ar ¨ar det praktiskt att definiera koefficienterna i ortogonalutvecklingen som en funktion av k,

c:=k->skal(P(k,x),f)/skal(P(k,x),P(k,x));

Sedan kan delsumman q_n av ortogonalutvecklingen skrivas

(10)

q:=n->sum(c(k)*P(k,x),k=0..n);

Plotta sen f och delsummor q_n samma diagram med t ex plot([f,q(5)],x=-1..1);

för n˚agra olika värden p˚a n, t ex n = 5, 7, 11, 13. Titta p˚a storleken av koefficienterna för n˚agra värden p˚a n, t ex evalf(c(10));, evalf(c(11));, evalf(c(20));, evalf(c(21));.

Gör samma med sin(10x) utbytt mot stegfunktionen f:=Heaviside(x);. Tag n = 5, 11, 25. Observera att liksom för Fourierserier uppkommer här ett Gibbs-fenomen. Titta p˚a storleken av koefficienterna med evalf.

Fr˚aga:Kan man g¨ora n˚agra iaktagelser om koefficienternas storlek i de b˚ada fallen?

5. F¨orberedelse: Skriv upp definitionen p˚a egenv¨arde och egenfunktion till en differentialoperator A.

Kontrollera att Legendrepolynomen ¨ar egenfunktioner till Legendres differentialoperator

− d

dx(1 − x²) d

dx. Skriv

d:=n->-diff((1-x^2)*diff(P(n,x),x),x)/P(n,x);

Skriv sen exempelvis

d(5);

och förenkla därefter med simplify(%); Pröva med n˚agra olika värden p˚a n, t ex 6, 7, 8, 10, 20.

Fr˚aga: Vad har detta med egenvärden att göra? Vad är egenvärdena uttryckta i n?

Stämmer detta med teorin? (Jämför sats S.4 sidan 340 i läroboken.)

4 Besselfunktioner

6. Besselfunktionerna J_νoch Y_νfinns naturligtvis i Maple, d¨ar de anropas med BesselJ(n,k);

resp BesselY(n,k);, se ?Bessel. ¨Aven nollst¨allena till Besselfunktionena finns i Maple.

De f˚as med kommandot BesselJZeros(n,k);. F¨or att slippa skriva s˚a mycket kan det vara praktiskt att inf¨ora kortare beteckningar, t ex

alias(J=BesselJ, Y=BesselY, alpha=BesselJZeros);

Rita n˚agra grafer, tex med plot([J(n,x)$n=0..5],x=0..15);. Man kan ocks˚a uppfatta J_ν(x) som en funktion av ν och x och rita en tredimensionell bild med

plot3d(J(v,x),x=0..10,v=0..10);

och vrida och vända p˚a denna. Lägg in lämpliga axlar.

Med kommandot matrix(10,9,(k,n)->evalf(alpha(n-1,k))); genereras en mas- sa nollställen. Jämför med tabellen i formelsamlingen.

(11)

7. Pröva om Maple klarar övning S.11. (Här behöver man använda simplify.)

8. Maple kan m˚anga formler för Besselfunktioner. Pröva att beräkna J_ν⁰(x). Pröva ocks˚a att uttrycka J1/2(x), J3/2(x), Y1/2(x) med elementära funktioner. Tänk p˚a att skriva J(3/2,x), inte J(1.5,x).

9. Extrauppgift. L¨os Bessels differentialekvation u⁰⁰(x) + 1

xu⁰(x) + (λ − ν²

x²)u(x) = 0

med dsolve (jämför sidan 6 i denna handledning). Här m˚aste du själv tänka p˚a att undersöka fallen λ = 0 och λ = 0, ν = 0. (Jämför lemma S.1 (sid 333) och övn S.13.) Lös även differentialekvationen

u⁰⁰(x) + 2

xu⁰(x) + (λ − `(` + 1)

x² )u(x) = 0 .

(Jämför (20) sid 340 och övning S.18.) Pröva speciellt ` = 0. Glöm inte fallet λ = 0.

10. Förberedelse: Repetera exemel S.3, sidan 334 i boken. Enligt detta exempel är Bes- selfunktionerna J0(α0kr) egenfunktioner till den singulära Sturm-Liouvilleoperatorn

Au = −1 r(ru⁰)⁰

DA=u ∈ C²[0, 1] | u begränsad nära 0, u(1) = 0 . Ange en skalärprodukt i vilken funktionernaJ⁰(α⁰_kr) är ortogonala.

Skriv in rätt viktsfunktion och intervallgränser. (Tänk p˚a att variabeln heter x i kommandot skal.) Kontrollera ortogonaliteten mellan J0(α0kr) för n˚agra olika k-värden med hjälp av skal. Plotta ocks˚a produkten J⁰(α⁰_kr)J⁰(α⁰mr)w(r), där w är viktsfunktionen,

¨over intervallet [0, 1].

Gör motsvarande undersökning för n˚agra Besselfunktioner J_n med n > 0. Slutsats?

11. Titta p˚a tredimensionella bilder av egenfunktionerna J_n(α_nkr) cos(nθ) till Laplaceo- peratorn med Dirichletvillkor i polära koordinater (exempel S.4 sid 335–337) för n˚agra olika värden p˚a n och k. För n = 2, k = 1 kan man skriva

plot3d([r*cos(v),r*sin(v),J(2,alpha(2,1)*r)*cos(2*v)],r=0..1,v=0..2*Pi);

Fast inspända cirkulära membran (se exempel 3.14 sid 108) har svängningsmoder av typ

u_nk = Jn(α_nkr) cos nθ sin cα_nkt .

S˚adana svängningsmoder kan man se röra sig i Maple med hjälp av kommandot animate3d.

Detta f˚ar man tillg˚ang till om man skriver

with(plots);

F¨or att se membranet sv¨anga i fallet n = 2, k = 1, och med c vald s˚a att cα_nk = 1, kan man skriva

(12)

animate3d([r*cos(v),r*sin(v),J(2,alpha(2,1)*r)*cos(2*v)*sin(t)], r=0..1,v=0..2*Pi,t=0..2*Pi);

Filmen startas och styrs fr˚an menyraden (starta genom att klicka p˚a rutan m¨arkt I).

Genom att klicka p˚a lämplig knapp kan man f˚a den att g˚a oavbrutet. Titta p˚a sväng- ningsmoder för n˚agra andra värden p˚a n och k. Man kan f˚a en jämnare rörelse genom att öka antalet indeningspunkter i tids- eller rumsled, se ?animate.

Vill man titta p˚a ett membran som inte är fast inspänt, utan kan röra sig fritt p˚a randen, p˚a s˚a sätt att det uppfyller ett homogent Neumannvillkor, s˚a kan man byta α_nk mot α⁰_nk (nollställen till J_n⁰). Dessa kan hämtas tabellen i formelsamlingen eller f˚as med fsolve. Pröva t ex

fsolve(diff(J(1,x),x)=0,x=1..3);

Här har vi angivit i vilket intervall det sökta nollstället ligger. För att kunna göra detta

¨ar det l¨ampligt att plotta funktionen J1⁰.

12. Extrauppgift:Utveckla funktionen f (x) = x, 0 < x < 1 i den ortogonala basen J0(α0kx).

Plotta delsummor S_N och funktionen f i samma diagram för n˚agra olika värden p˚a N . Gör sedan samma sak med f (x) = 1 − x. Märks n˚agon skillnad? Vad kan det bero p˚a?

F¨orberedelser:

Skriv upp formeln f¨or skal¨arprodukt i L2([0, L]): (f |g) =

Antag att ortogonalutveckligen f =P∞

k=1c_kϕ_k gäller. Beskriv Fourierkoefficienterna c_k med hjälp av skalärprodukt: c_k =

H.22Laguerrepolynomen Ln(x) är ortogonala i L²(I, w), där I = [0, ∞) och w(x) = e⁻^x. De tre första är

L⁰(x) = 1, L¹(x) = 1 − x, L²(x) = x²− 4x + 2.

Bestäm med hjälp härav det polynom p = p_α av andra graden som minimerar integralen Z ^∞

0

(e⁻^αx− p(x))²e⁻^xdx, d¨ar konstanten α > −1/2.

Skriv upp en anv¨andbar skal¨arprodukt: (f |g) =

Skriv med projektionsformeln upp en formel f¨or p: