# Dubbelintegraler med Maple
# 
# Påminnelse om enkelintgeraler: Integralen av funktionen t->f(t) från
# t=a till t=b beräknas med kommandot int. Ex:
> int(t^2+sin(t),t=a..b);
> 
# Om vi istället skriver Int (med stort I) får vi integralen skriven med
# integraltecken:
> Int(t^2+sin(t),t=a..b);
> 
# Detta är ett bra sätt att kontrollera att det är rätt integral man
# beräknar.
# 
# Som framgår av boken återförs beräkning av dubbelintegraler på
# upprepade enkelintegraler. I Maple utförs dessa (naturligt nog) genom
# upprepade användningar av kommandot int. Liksom förut kan vi använda
# Int för att få integralen med integraltecken. 
# 
# Exempel Med  f(x,y)=(x-y)*sin(y) och integrationsområdet givet av
# 0<x<2, 0<y<Pi (en rektangel)  har vi
> Int(Int((x-y)*sin(y),y=0..Pi),x=0..2);#   
> Int(int((x-y)*sin(y),y=0..Pi),x=0..2);
> int(Int((x-y)*sin(y),y=0..Pi),x=0..2);
> int(int((x-y)*sin(y),y=0..Pi),x=0..2);
# Observera att kombinationen int (Int) ger samma resultat som Int(Int).
> 
> 
# Exempel 4.6 Vi börjar med att rita en figur över integrationsområdet.
> with(plots):
> implicitplot({x=1,x=3,y=-2*x/3+4,y=x/3},x=0..5,y=0..5);
# Vi kan nu beräkna integralen genom att integrera först i y- och sedan
# i x-led.
> I=Int(Int(3*x+18*y,y=x/3..-2*x/3+4),x=1..3);
> I=int(int(3*x+18*y,y=x/3..-2*x/3+4),x=1..3);
# Anmärkning: Att först integrera i x-led är mindre lämpligt eftersom vi
# då får tre olika uttryck för den övre integrationsgränsen x_max(y) 
# beroende
# på vilket intervall y tillhör.
# 
# Exempel på en trippelintegral
# 
# Uppgift 4.33. 
# Eftersom gränserna i x och y är oberoende av z verkar det rimligt att
# börja integrera i z-led.
> Int(Int(Int(x+y,z=0..exp(x*y)),y=0..1),x=0..1);
> Int(Int(int(x+y,z=0..exp(x*y)),y=0..1),x=0..1);
# Med Maple är det inget problem att integrera hela uttrycket först med
# avseende på y och sedan m a p x (eller tvärtom--
# x och y spelar uppenbarligen samma roll).
# 
> Int(int(x*exp(x*y)+y*exp(x*y),y = 0 .. 1),x = 0 .. 1);
> int((exp(x)*x^2-exp(x)+exp(x)*x-x^2+1)/x^2,x = 0 .. 1);
# Observera dock att den sista integralen skulle vara mycket
# ansträngande att beräkna för hand. (Det tar ett tag att direkt från
# funktionsuttrycket
# se att integranden faktiskt är en kontinuerlig och begränsad
# funktion--vilket är klart för oss eftersom den är resultatet av att
# integrera en kontinuerlig funktion av x och y över ett ändligt
# y-intervall)  En primitiv funktion som försvinner i x=0 ges av 
# 
> int((exp(x)*x^2-exp(x)+exp(x)*x-x^2+1)/x^2,x=0..X);
> 
# Dessa svårigheter kan undvikas om vi observerar att dubbelintegralen
# är en summa av
# 
> I1:=Int(Int(x*exp(x*y),y=0..1),x=0..1);
# och
> I2:=Int(Int(y*exp(x*y),y=0..1),x=0..1);
# där I1 och I2 enbart skiljer sig åt genom att variablerna x och y bytt
# roll (observera att detta inte ändrar integrationsområdets utseende). 
# Det gäller alltså I1=I2. Men den inre integralen i I1 har en mycket
# enkel form
> I1:=Int(int(x*exp(x*y),y=0..1),x=0..1);
> I1:=int(int(x*exp(x*y),y=0..1),x=0..1);
#                     
#                                            Övergång till polära
# koordinater
# 
# Detta kommer att diskuteras mer senare.
# 
# Antag att vi vill beräkna integralen av f(x,y)=x+y över
# halvcirkelskivan x^2+y^2<4, x>0. 
> implicitplot({x=0,x^2+y^2=4},x=0..3,y=-3..3,scaling=constrained);
# 
# 
# 
# Upprepad integration är här inget problem för Maple så snart man
# observerat att gränserna i y-led fås genom att man löser ut y ur
# x^2+y^2=4,
# dvs y_min=-sqrt(4-x^2) och y_max=sqrt(4-x^2):
# 
> SöktaIntegralen:=Int(Int(x+y,y=-sqrt(4-x^2)..sqrt(4-x^2)),x=0..2);
# 
> SöktaIntegralen:=Int(int(x+y,y=-sqrt(4-x^2)..sqrt(4-x^2)),x=0..2);
> SöktaIntegralen:=int(int(x+y,y=-sqrt(4-x^2)..sqrt(4-x^2)),x=0..2);
# Vid handräkning (och för mer allmänna cirkelsektorer) är det dock
# fördelaktigt att gå över till polära koordinater x=r*cos(v),
# y=r*sin(v), se Sats 5B på sidan 287. Härvid övergår dxdy i
# dubbelintegralen i rdrdv. (Observera faktorn r  som kommer från
# funktionaldeterminanten för variabelbytet!) 
# I polära koordinater blir halvcirkelskivan rektangel 0<r<2,
# -Pi/2<v<Pi/2 och integranden blir (r*cos(v)+r*sin(v))*r. (Det sista
# r:et kommer
# från funktionaldeterminanten.
# 
# Vi får 
> 
#   
> SöktaIntegralen:=Int(Int((r*cos(v)+r*sin(v))*r,r=0..2),v=-Pi/2..Pi/2);
>  SöktaIntegralen:=Int(int((r*cos(v)+r*sin(v))*r,r=0..2),v=-Pi/2..Pi/2)
> ;
>  SöktaIntegralen:=int(int((r*cos(v)+r*sin(v))*r,r=0..2),v=-Pi/2..Pi/2)
> ;
>