Transcript Rapport TMA682 110914.pdf

Projekt
Finit Element-l¨osare
Emil Johansson, Simon Pedersen, Janni Sund´en
29 september 2011
Chalmers Tekniska H¨ogskola
Institutionen f¨or Matematik
TMA682 Till¨ampad Matematik
1
Inledning
M˚
anga naturliga fenomen kan modelleras med differentialekvationer. Ofta ¨ar dessa mycket sv˚
ara eller om¨
ojliga att l¨
osa analytiskt, varf¨or olika approximativa l¨osningsmetoder
anv¨ands. En av dessa ¨
ar Finita element-metoden (FEM). Denna rapport behandlar ett
projekt best˚
aende av kodning av en FEM-l¨osare i MatLab samt anv¨andandet av denna p˚
a tv˚
a v¨
alk¨
anda randv¨
ardesproblem, diffusions-reaktionsekvationen och diffusionskonvektionsekvationen. Den erh˚
allna approximativa l¨osningen j¨amf¨ors med en analytisk
l¨osning av problemet.
2
Teori
2.1
Finita element-metoden
Differentialekvationens l¨
osning kan approximeras med en styckvis linj¨ar, kontinuerlig,
funktion. Detta kallas kontinuerlig Galerkin av grad 1, cG(1). D˚
a denna approximativa l¨osning ¨
ar (styckvis) linj¨
ar kan den beskrivas som en linj¨arkombination av styckvis
linj¨ara basfunktioner, s˚
a kallade hattfunktioner. Att endast beh¨ova unders¨oka dessa basfunktioner ¨
ar grunden i l¨
osningsmetoden. Som exempel har vi h¨ar behandlat diffusionsreaktionsekvationen, men proceduren ¨ar liknande f¨or andra randv¨ardesproblem.
(
− u00 (x) + u(x) = 1, 0 < x < π
(DE)
(1)
u0 (0) = 1, u(π) = 0
F¨or att problemet ska kunna behandlas, konverterar vi differentialekvationen till motsvarande variationsformulering genom att multiplicera (DE) med en testfunktion v(x) och
integrera ¨
over det relevanta intervallet. F¨or v˚
art exempelproblem (ekv. 1) blir denna:
Finn en funktion u(x) ∈ V s˚
adan att:
Z π
Z π
(−u00 (x) + u(x))v(x) dx =
(1)v(x) dx
0
0
Z π
1
∀ v(x) ∈ V := H (0,π) := {w :
w(x)2 + w0 (x)2 dx < ∞, w(π) = 0}
(VF) (2)
0
Via partiell integration erh˚
aller vi f¨oljande f¨orenklade formel:
x=π Z
− u0 (x)v(x) + u(x)v(x)
+
π
u0 (x)v 0 (x) dx +
0
x=0
Z
0
π
u(x)v 0 (x) dx =
Z
π
v(x) dx
0
(3)
∀ v(x) ∈ V
Ins¨attning av randdata l¨
amnar endast kvar integralerna, d˚
a v(pi) = 0 och u0 (0) = 0.
Z π
Z π
Z π
0
0
0
u (x)v (x) dx +
u(x)v (x) dx =
v(x) dx, ∀ v(x) ∈ V
(4)
0
0
0
Testrummet V ¨
ar allts˚
a alla funktioner som ¨ar begr¨ansade och har begr¨ansade derivator, samt uppfyller randvillkoren f¨
or (DE). Detta testrum ¨ar o¨andligt dimensionellt och
testfunktionerna ¨
ar o¨
andliga till antalet. F¨or att hitta en approximativ l¨osning beh¨over
vi inskr¨anka oss till ett ¨
andligt dimensionellt rum. Detta ¨ar sj¨alva finita element-metoden.
1
L˚
at Th : 0 = x0 < x1 < · · · < xi < xi+1 < · · · < xN = π vara en likformig partition av
det relevanta intervallet med delintervall¨angden h.
Finn en approximativ l¨
osning U (x) ∈ Vh s˚
adan att:
Z π
Z π
Z π
0
0
0
v(x) dx,
U (x)v (x) dx =
U (x)v (x) dx +
(5)
0
0
0
∀ v(x) ∈ Vh := {w : w ¨
ar styckvis linj¨ar p˚
a Th och kontinuerlig, w(π) = 0}
Den approximativa l¨
osningen sammanbinder nodpunkterna (x1 , xi osv) med styckvis
linj¨ara segment. Vi till˚
ater h¨
armed v˚
ar funktion att avvika fr˚
an den analytiska l¨osningen
s˚
a l¨ange den ¨
ar lika i nodpunkterna.
Eftersom U (x) ¨
ar styckvis linj¨
ar kan den skrivas som en summa av styckvis linj¨ara
N
basfunktioner {ϕj (x)}j=1 med koefficienterna ξ0 , . . . , ξN .
U (x) =
N
X
ξj ϕj (x) =⇒ U 0 (x) =
j=1
N
X
ξj ϕ0j (x)
(6)
j=1
Ins¨attning i (5) ger oss f¨
oljande:
N Z
X
j=1
π
ξj ϕ0j (x)v 0 (x) dx +
π
Z
ξj ϕj (x)v 0 (x) dx =
π
v(x) dx,
∀ v(x) ∈ Vh
(7)
0
0
0
Z
Vi v¨
aljer nu testfunktionerna som hattfunktionerna ϕ1 (x), . . . ,ϕ(N −1) (x). Eftersom
ingen homogen dirichletdata finns i intervallets startpunkt kr¨avs ocks˚
a en halv hattfunktion som bas i x0 . V˚
art problem ¨
ar nu formulerat som f¨oljande:
N Z
X
j=1
π
ϕ0j (x)ϕ0i (x) dx
ξj +
0
Z
π
ϕj (x)ϕ0i (x) dx
0
Z
ξj =
π
ϕi (x) dx,
i = 1, 2, . . . , N
0
(8)
Ekvationen kan nu skrivas om p˚
a matrisform. Vi f˚
ar f¨oljande samband:
Z π
N
S = sij i,j=1 , sij =
ϕ0j (x)ϕ0i (x) dx
0
Z π
N
M = mij i,j=1 , mij =
ϕj (x)ϕ0i (x) dx
0
 
 
ξ1
b1
Z π
 b2 
 ξ2 
 
 
ϕi (x) dx
ξ= . 
b =  .  bi =
 .. 
 .. 
0
ξN
(9)
bN
Om differentialekvationen hade inneh˚
allit en termRav f¨orsta ordningen, hade en konπ
vektionsmatris uppst˚
att vars element a¨r beroende av 0 ϕj (x)ϕi (x) dx.
Med hj¨
alp av resonemang kring hattfunktionernas st¨od och ¨overlapp kommer vi fram
till sambandet att
2



sij






s
ij




sij





sij



mij






m
ij




mij





mij
1
, i=j=1
h
2
= , i=j>1
h
−1
=
, ki − j| = 1
h
= 0, ki − j| > 1
=
2h
, i=j=1
6
4h
=
, i=j>1
6
1
= , ki − j| = 1
h
= 0, ki − j| > 1
=
Det vill s¨
aga, matriserna S och M f˚
ar f¨oljande utseende:



1 −1 0 · · · 0
2 1



.
.
−1 2 −1 0
1 4
. 

1
h


.
.
..
.
.
..
..
..
..  M = 
S= 0

.

h
6 0
 ..
 ..

..
 . ···
. · · ·
. 2 −1
0 · · · 0 −1 2
0 ···
(10)
0
···
1
..
.
..
.
0
..
.
0
4
1

0
.. 
.

.. 
.


1
4
(11)
Eventuell konvektionsmatris ¨
ar antisymmetrisk och har nollor p˚
a huvuddiagonalen.
Den ser ut som f¨
oljande (notera avsaknandet av h-beroende):


0
1
0 ··· 0

.
−1 0
1
0 .. 


1

.
.
.
.
.
.
.
.
(12)
K= 0
. .
.
.

2
 ..

..
 . ···
. 0 1
0 · · · 0 −1 0
Styvhets- och massmatriserna ¨
ar i allt v¨asentligt lika f¨or alla problem, endast avvikande i elementen (1,1) och (N,N ) om randvillkoret i f¨orsta respektive sista punkten ¨ar
nollskilda. D˚
a konvektionsmatrisens huvuddiagonal endast inneh˚
aller nollor p˚
averkas den
inte av eventuella halvhattfunktioner. De tre matriserna skalas ocks˚
a efter l¨angden p˚
a
00
0
varje delintervall (h) och koefficienterna framf¨or u (x), u (x) och u(x), vilket enkelt ses i
variationsformuleringen f¨
or problemet i fr˚
aga.
Med liknande resonemang f˚
as elementen i lastvektorn b.
R π Arean under varje hattfunktion
ojden 1 och basen 2h. Integralen 0 ϕi (x) dx motsvarar denna area
¨ar en triangel med h¨
och evalueras d¨
arf¨
or till h, utom i x0 d¨ar vi har en halv hattfunktion, varf¨or basen blir h
och integralens v¨
arde d¨
arf¨
or halveras. Lastvektorn ser allts˚
a ut som f¨oljande:
1
2
1
 
 
b = h 1
 .. 
.
1
(13)
Att finna den approximativa l¨
osningen till v˚
art randv¨ardesproblem inneb¨ar nu att l¨osa
matrisekvationen (S + M )ξ = b.
2.2
Finit element-l¨
osare i MatLab
MatLab-koden grundar sig p˚
a den inbyggda funktionen f¨or l¨osning av matrisekvationer av
typen Ax = b. Styvhets-, mass- och konvektionsmatriserna byggs upp genom iteration
3
och korrigeras f¨
or halvhattfunktioner d¨ar randvillkoren a¨r nollskilda. Detsamma g¨aller
lastvektorn b. De olika matriserna summeras sedan, och ekvationen (S + M + K)ξ = b
l¨oses med MatLabs backslash-operator. Koden i sin helhet finns i bilaga A. Randv¨ardesproblemen som l¨
oses ¨
ar diffusions-reaktionsekvationen
(
− u00 (x) + u(x) = 1, 0 < x < π
(14)
u0 (0) = 1, u(π) = 0
och diffusions-konvektionsekvationen

 − εu00 (x) + 1 u0 (x) = 1,
2

u(0) = u(π) = 0
0<x<π
(15)
d¨ar ε ¨ar en liten, konstant, godtycklig parameter, i koden satt som ε = 3.
Figur 1: Diffusions-reaktionsekvationens analytiska l¨osning samt cG(1)-l¨osning med antalet delintervall N = 4 respektive N = 13. Notera approximationens n¨arhet till den analytiska l¨osningen
redan vid l˚
aga v¨
arden p˚
a N. Approximationen ¨ar fixerad vid y = 0 i x = π f¨or att svara mot
dirichletdatan.
4
Figur 2: Diffusions-konvektionsekvationens analytiska l¨osning samt cG(1)-l¨osning med antalet
delintervall N = 4 respektive N = 13. Notera att den approximativa l¨osningen h¨ar ¨ar fixerad
i b¨
agge ¨
andpunkterna, vilket svarar mot dirichlet-randvillkoren.
A
A.1
Bilaga: MatLab-kod
Diffusions-reaktionsekvationen
%% D i f f u s i o n s r e a k t i o n s e k v a t i o n e n
clc ;
clf ;
a
c
N
h
x
=
=
=
=
=
0; % Intervallet :
pi ; % 0 < x < p i
3 ; % An tal d e l i n t e r v a l l = N+1 ( D e t t a g e r b ¨
a t t r e approximation )
( c−a ) / (N ) ; % D e l i n t e r v a l l ¨
angden
linspace ( a , c , N+1); % Nodpunkter
% Analytisk l¨
osning
u = @( x ) (−exp ( pi−x)−exp ( pi+x)+1+exp ( 2 ∗ pi ) ) . / ( 1 + exp ( 2 ∗ pi ) ) ;
f p l o t ( u , [ a c ] , ’ k− ’ ) ;
hold on ;
b = o n e s (N, 1 ) ; % L a s t v e k t o r n
% −− S t y v h e t s m a t r i s −−
S = zeros (N,N ) ;
for i = 1 :N
for j = 1 :N
i f ( i==j ) % Phi−prim ¨
overlappar h el t
S( i , j ) = 2;
e l s e i f ( abs ( i −j ) == 1 ) % H a l v t ¨
overlapp
S ( i , j ) = −1;
end
5
end
end
S ( 1 , 1 ) = 1 ; % H a l v h a t t p˚
a x0−x1
S=1/h∗S ;
% Massmatris
M = zeros (N,N ) ;
for i = 1 :N
for j = 1 :N
i f ( i==j ) % Som ovan f a s t Phi
M( i , j ) = 4 ;
e l s e i f ( abs ( i −j ) == 1 ) % . . .
M( i , j ) = 1 ;
end
end
end
M( 1 , 1 ) = 2 ; % H a l v h a t t p˚
a x0−x1
M=h/6∗M;
b(1) = 1/2; % Halvhatt !
b = h∗b ;
x i = ( S+M) \ b ; % L¨
o s A∗ x i = b
U = zeros (N+1, 1 ) ; % K o r r i g e r a f ¨
o r r a n d d a t a U( p i ) = 0
U( 1 :N) = x i ;
plot ( x , U, ’ r .− ’ ) ;
xlabel ( ’ x ’ ) ;
ylabel ( ’ y ’ ) ;
legend ( ’ y = u ( x ) ’ , ’ y = U( x ) ’ ) ;
A.2
Diffusions-konvektionsekvationen
%% D i f f u s i o n s −k o n v e k t i o n s e k v a t i o n e n
clc ;
clf ;
a = 0; % Intervallet :
c = pi ; % 0 < x < p i
N = 3 ; % An tal d e l i n t e r v a l l = N+1 ( D e t t a g e r b a
¨ t t r e approximation )
h = ( c−a ) / (N+1); % D e l i n t e r v a l l ¨
angden
x = linspace ( a , c , N+2); % Nodpunkter
e = 3; % Godtyckligt , l i t e t , epsilon
b = o n e s (N, 1 ) ;
% Analytisk l¨
osning
u = @( x ) 2 ∗ ( ( exp ( pi / ( 2 ∗ e )) −1)∗ x−pi ∗exp ( x / ( 2 ∗ e ))+ pi ) / ( exp ( pi / ( 2 ∗ e ) ) − 1 ) ;
f p l o t ( u , [ a c ] , ’ k− ’ ) ;
6
hold on ;
% Styvhetsmatris
S = zeros (N,N ) ;
for i = 1 :N
for j = 1 :N
i f ( i==j ) % Som ovan
S( i , j ) = 2;
e l s e i f ( abs ( i −j ) == 1 ) % Som ovan
S ( i , j ) = −1;
end
end
end
S=e /h∗S ; % Ta h¨
a nsyn t i l l e p s i l o n
% Konvektionsmatris
C = zeros (N,N ) ;
for i = 1 :N
for j = 1 :N
i f ( ( i −j ) == 1 ) % Under d i a g o n a l e n ( rad > k o l o n n )
C( i , j ) = −1;
¨
e l s e i f ( ( j −i ) == 1 ) % Over
d i a g o n a l e n ( rad < k o l o n n )
C( i , j ) = 1 ;
end
end
end
C=1/2∗1/2∗C ; % Gl¨
om i n t e k o e f f i c i e n t e n f r a m f o
¨r u ’ ( x ) !
b ( 1 ) = 1 ; % Inga h a l v h a t t a r
b (N) = 1 ; % . . .
b = h∗b ;
x i = ( S+C) \ b ;
U = zeros (N+2, 1 ) ; % K o r r i g e r a f o
¨ r r a n d d a t a U( 0 ) = U( p i ) = 0
U( 2 :N+1) = x i ;
plot ( [ x ] , U, ’ r .− ’ ) ;
xlabel ( ’ x ’ ) ;
ylabel ( ’ y ’ ) ;
legend ( ’ y = u ( x ) ’ , ’ y = U( x ) ’ ) ;
7