Metody numeryczne

dr hab. Piotr Fronczak

Metody numeryczne

Wykład nr 10

Równania różniczkowe cząstkowe

Równania różniczkowe cząstkowe (RRC)

• liczba zmiennych 2

 

 







, ,

, ,

0 ,

, , , , , , ,

, ,

2 2

2

y x u u

x u u

y u u

x u u

u u u u u u y x F

y x u u

xy xx

y x

yy xy xx y x





 



 



 



 





• rząd równania: rząd najwyższej pochodnej

0 0

0

3 3













yyy x

y xx

y x

u u

u u

u b u

• charakterystyka: liniowe, quasi-liniowe, nieliniowe

     

   

   

   

 

           

     

   

   

  

x y xx xy yy



y x

y x

yy xy

xx y x

y x

y x

y x

u u u u u u y x

u u u y x

y x

C B A

F u

E u

D u

C u

B u

A

u u

x u

u u

bu u

u u u y x

u y x

y x

c b a

, , , , , , , :

, , , , :

, :

: ： , : , :

0 :

: :

0

0 0

, , , ,

, , , , ,

2

2



























































Nieliniowe liniowe Quasi

Liniowe ...

...

...



 

^{ u  b}

 

^{ u  c}

 

^{  0}

a _{x y}

Nieliniowe

liniowe Quasi

Liniowe ...

...

...



     

   

   

   

 

⁰

0 0

, , , ,

, , , , ,

2

2



































y x

y x

y x

y x

u u

x u u u

bu u

u u u y x

u y x

y x

c b a

• Skupmy się na RRC co najwyżej drugiego rzędu – najpopularniejsze w fizyce

• Ogólna postać RRC drugiego rzędu

B²-4AC Kategoria Przykład

< 0 eliptyczne Równanie Laplace’a

= 0 paraboliczne Równanie przewodnictwa ciepła

>0 hiperboliczne Równanie falowe









2 2

2

2 0

T x

T

 y 

k T x

T t









2 2 









2

2 2

2 2

1 y

x c

y

 t

) , , , , ( )

, ( )

, ( )

,

( ₂

2 2

2 2

y y x

u x u u y G

y u x y C

x y u x x B

y u x

A 





 



 





 





) , , , ,

( u u u x y G

Cu Bu

Au

_xx



_xy



_yy



_{x y}

Motywacja dla takiej klasyfikacji

eliptyczne C

y u ne x

parabolicz C

x u zne

hiperbolic C

y u

x  



 







 



 



 







 



 







 





2 2 2

2 2

2 2

2 2

2

;

; Najprostsze rozwiązania:

Również dla bardziej skomplikowanych równań lokalne własności rozwiązanie zależą od znaku wyrażenia B²-4AC.



^{x y}



^{hiperbola u C x parabola u C}



^{x y}



^elipsa

u C ^{2 2 2 2 2}

; 4

; 2

4    



Zagadka: Skategoryzuj równanie Schrodingera.

 0



_y

xx

iCu

u

• Zakładamy h_x = h_y = h [siatka kwadratowa]

• u(x_i,y_j) = u_i,j

• u(x_i+h, y_j+h) = u_i+1,j+1,



i j



x

 u

i j

u

i j

u

i j



y h x

u 1

_{2 1,}

2

_{, 1,}

,

'' 

_

 

_

x

_i-1

x

_i

x

_i+1

y

_j

y

_j-1

y

_j+1

u

_i,j ^hx

h_y

Eliptyczne RRC - dwuwymiarowe zagadnienie brzegowe

2

0

2 2

2





 





y u x

Równanie Laplace’a

u

  1

2



, 1

2

, , 1



,

'' 

_{i j}



_{i j}



_{i j}

j y

i

u u u

y h x u

   

   

 

 ⁴  ⁰

1

2 1 2

, '' ,

''

1 , ,

1 ,

1 , ,

2 1

1 , ,

1 , ,

1 ,

, 2 1











































j i j

i j

i j

i j

i

j i j

i j

i j

i j

i j

i

j y x i

j i

u u

u u

h u

u u

u u

u h u

y x u y

x

u

Sumując wyrazy wyznaczamy u_i,j: _,



_{1, , 1 1, , 1}



4 1









  



_{i j i j i j i j}

j

i

u u u u

u

1 2 3 4

u = 0

u = 0 u = 0

u = 1

Przykład: siatka 2x2

) 1

0 (

) 0 1 (

) 0

0 (

) 0 0

(

3 2

41 1 4

4 4 1

1 1 3

1 4 4

1 1 2

3 4 2

1 1 1

m m

m

m m

m

m m

m

m m

m

u u

u

u u

u

u u

u

u u

u









































Przyjmujemy początkowe przybliżenie u₁¹ u¹₂  u₃¹ u¹₄  ₀

Korzystamy z metody iteracyjnej Jacobiego (lub Gaussa-Seidla – szybsza zbieżność)

Korzystając z metod dokładnych (np. dekompozycja LU) musimy ułożyć macierz o rozmiarze liniowym n x n.

Wniosek: Musimy wprowadzić indeksowanie równań odpowiadających punktom dwuwymiarowej siatki n x n:

) ,

2

(

2 2

2

y x y f

u x

u  



 



Przykład: równanie Poissona



 u

i j

u

i j

u

i j

u

i j

u

i j

 f

i j

h 1

²



_^1,



^, _¹

 4

^,



_^1,



^, _¹



^,

...,

2

, 2 ,

1 n

P P 

j

i , )  (

(na przykład wierszami)

i n j

P  (  1 )  

 ^u

_P

^u

_{P n}

^u

_P

^u

_P

^u

_{P n}

 ^f

_P

h 1

₂



_1



_

 4 

_1



_





0

  



potencjał

gęstość ładunku

Przykład:

Powierzchnia potencjału przy losowo rozmieszczonej gęstości ładunku

 ^u

_P

^u

_{P n}

^u

_P

^u

_P

^u

_{P n}

 ^f

_P

h 1

₂



_1



_

 4 

_1



_



Paraboliczne i hiperboliczne RRC

Bezwymiarowe równanie przewodnictwa ciepła

t

xx

u

u

t t x t u

x u





 

 ( , )

) ,

2

(

z warunkiem początkowym (dla t = 0)

L x

x g x

u ( , 0 )  ( ) 0  

gdzie L = 1 – szerokość dziedziny rozwiązań.

0 )

( )

, (

0 )

( )

, 0 (









t t

b t

L u

t t

a t

u

gdzie a(t) i b(t) są funkcjami jedynie czasu, a g(x) zależy jedynie od położenia x.

Warunki brzegowe:

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

Dyskretyzacja przestrzeni rozwiązań

czas (indeks j)

przestrzeń

(indeks i) x

t

 t

 x

 



















1 0

0 n i

x i x

j t

j t

i

Położenie w węźle i w chwili j: j

Zamieniając pochodne na różnice otrzymujemy

t

t x u t

x u t

t x u t

t x

u

_t

u

^{i j i j}



 







 ( ,  ) ( , ) ( ,

^1

) ( , )

z błędem O(t),

2

1

1

, ) 2 ( , ) ( , )

(

x

t x u t

x u t

x

u

_xx

u

^{i j i j i j}







^



^

z błędem O(x²).

2

1 1

1

2 x

u u

u u

t u u u

j i j

i j

xx i

j i j

t i





 



 





Dalej upraszczając 

Zatem zdyskretyzowane równanie ma postać

2

1 1

1

2

x

u u

u t

u

u

_i^j _i^j _i^j _i^j _i^j





 





_{ }



Definiując r = t / x², otrzymujemy schemat jawny Eulera:

j i j

i j

i j

i

ru r u ru

u

^1



_1

 ( 1  2 ) 

_1

- punkt uwzględniany przy obliczaniu różnicy czasowej

- punkt uwzględniany przy obliczaniu różnicy przestrzennej

// mamy n+1 węzłów przestrzennych siatki // warunek brzegowy

u[0] = u[n] = 0.0;

// warunek początkowy for (int i=1; i<n; i++) {

x = i*dx;

u[i] = func(x);

};

// pętla czasowa

for (int t=1; t<=tsteps; t++) for (int i=1; i<n; i++)

unew[i] = r*u[i-1] + (1-2*r)*u[i] + r*u[i+1];

u  unew;

Algorytm schematu jawnego

Problem: warunek stabilności schematu

2 1

2





 x

t

Czyli jeśli podzielimy domenę [0,1] na sto podprzedziałów (x = 0.01), to t ≤ 5·10^-5.

j i j

i j

i j

i

ru r u ru

u

^1



_1

 ( 1  2 ) 

_1

Analiza stabilności

Rozważmy schemat jawny dla równania ciepła:

  ¹

2



1

²

1

 ^{( 1 )}

1 2

2

n j n

j n

j n

j n

j

u u u

t x u u

x u t

u









 

 

 



 







 





Zatem nasz schemat ma postać:

 

2



1

²

1

 ^{( 2 )}

1 n

j n

j n

j n

j n

j

u u u

x u t

u

^ _

 

_



 



Niech

D

będzie dokładnym rozwiązaniem równania (1).

Niech

N

będzie numerycznym rozwiązaniem równania (1).

Zatem błąd schematu

D

 N 



Napiszmy równanie stabilności metody numerycznej, które opisze ewolucję błędu w trakcie obliczania kolejnych kroków czasowych.

stabilność nierośnie

ość niestabiln

rośnie



 

Równanie stabilności będziemy badać za pomocą analizy Fouriera (metoda von Neumanna).

Rozwiązanie numeryczne możemy zapisać jako

N = D + ε

(3)

Podstawiając (3) do (1), otrzymujemy

 

^{2 1 1}

1 1

1

1

2 2

x

D D

D t

D

D

ⁿ_j ⁿ_j ⁿ_j ⁿ_{j j}^{n n}_{j j}^{n n}_{j j}^{n n}_j











 









^ _{   }



    

    _ ^ _



 









 

 





 









 



 





^ _{   }



2

1 1

2

1 1

1

1

2 2

x x

D D

D t

t D

D

ⁿ_j ⁿ_j



ⁿ_j



ⁿ_j ⁿ_j ⁿ_j ⁿ_j



ⁿ_j



ⁿ_j



ⁿ_j

Ponieważ z równania (1)

  _ ^ _



 









 





_{ }



2

1 1

1

2

x

D D

D t

D

D

ⁿ_j ⁿ_j ⁿ_j ⁿ_j ⁿ_j

  ^{( 4 )}

2

2

1 1

1

 





 









 





_{ }



t x

n j n

j n

j n

j n

j

   

Zatem



ε (x,t)

x

Rozważmy rozkład błędu w pewnym kroku czasowym.

Błąd

ε(x,t)

można zapisać w postaci szeregu Fouriera:

 

) :

( , sin

cos

) 5 ( )

( ,

falowa liczba

m m

m x

ik

m

x ik m

k x

k i

x k e

e t b t

x

m

m





 



Czyli równanie na błąd jest takie same jak równanie na funkcję

u

^.

  ^{x , t  b}

_m

^{( t ) e}

^ikmx

^{( 6 )}



Ponieważ równanie różnicowe na błąd



jest liniowe, zachowanie każdego wyrazu szeregu jest podobne do zachowania całego szeregu. Zatem wystarczy rozważyć wzrost błędu typowego wyrazu

Czasową zależność błędu uwzględnimy, przyjmując, że amplituda błędu

b

_m jest funkcją czasu. Ponieważ błąd rośnie lub maleje zwykle wykładniczo z czasem, możemy zapisać

  ^{x  , t e}

^at

^e

^ikmx

^{( 7 )}



gdzie

k

_m jest rzeczywiste, ale

a

może być zespolone.



^{t t}



^{ik x at ik x}



^{at ik}



^{x x}

 ²

^{at ik x at ik}



^{x x}

  ^{( 8 )}

a

e

_m

e e

_m

r e e

_m

e e

_m

e e

_m

e

^

 

^

 

^

 

x ² r t

gdzie



 

Podstawiając (7) do (4), otrzymujemy

Dzieląc (8) obustronnie przez

e

^at

e

^ikmx



^{t t}



^{ik x at ik x}



^{at ik}



^{x x}

 ²

^{at ik x at ik}



^{x x}

  ^{( 8 )}

a

e

_m

e e

_m

r e e

_m

e e

_m

e e

_m

e

^

 

^

 

^

 ²  ^{( 9 )}

1

^{ik x ik x}

t

a

r e

_m

e

_m

e

^

 

^

 

^{ }

cos 2



  ^e

ⁱ

^{ e}

^i

Korzystając z zależności

) 1 (cos

2

1  



 r 

e

^{a t}

otrzymujemy

gdzie

x k

_m



 

2 cos 1

sin

²

 2   

Korzystając z zależności

sin 2 4

1

²



r e

^at

 

otrzymujemy

n j n

G

j





^1



Definiujemy współczynnik wzmocnienia błędu

Oczekujemy, by błąd nie narastał z kroku na krok

1

|

| G 

Zatem

, 1

1



 n j n j





2 1 sin 4

1 

²

 

r

  ^{x  , t e}

^at

^e

^ikmx



t a n

j n

j

e

^



 



¹

2 1 sin 4

1 

²

 

r

2 1 sin 4 1

1  

²





_r



2 0 1

) 2 ( ) 1

(   

r



0 2 0

sin 4 2 1

sin 4 1 ),

1

( 

r ²

  

r ²

  

r



2 1 2

1 sin 2

sin 2 4 1 1 ),

2

(   

r ²

 

r ²

  

r



2 1

2





 x

t

Warunek zbieżności schematu jawnego (1)

n j n

G

j





^1



krok czasowy

wykładnik potęgi Przypomnijmy sobie, że współczynnik wzmocnienia błędu

n j n

j

G 



^1



0 1 1

j n n

j

G 



^



^

Czyli:

Zwróćmy uwagę, że nasz schemat jawny (4)

  

ⁿj



n j n

j n

j n

j n

j n

j n

j n

j r

t ¹ x ₂ ¹ _{1 1}

1

2 2











    



 









 



       





 

 

 





 

 

 







 

 

 





 

 

 







 

 

 





 

 

 





 

 

 





 

 

 



















n N n n

n N n n

n N n n

G r

r r r r

r r

r

r r



















2 1

1 1 2

1 1 2

1

2 1 2

1 2

1 2

1

możemy przedstawić w postaci równania macierzowego

Równanie własne z wartościami własnymi

G

Kolejne przybliżenia x⁽⁰⁾, x⁽¹⁾ , x⁽²⁾ utworzą ciąg wektorów. Jeżeli istnieje granica tego ciągu, wtedy jest ona rozwiązaniem układu równań liniowych.

Ciąg wektorów musi być zatem ciągiem zbieżnym.

Tw. Ciąg określony wzorem (

*

) przy dowolnym wektorze x⁽⁰⁾ jest zbieżny wtedy i tylko wtedy, gdy ()<1.

() – promień spektralny macierzy  = max |_i|, _i – wartości własne macierzy . Jako rozwiązanie początkowe, obiera się dowolny wektor (np. wektor

zerowy) i oblicza się kolejne iteracje:

Metoda Jacobiego

Przypomnienie z wykładu nr 3

Zatem mamy układ równań iterowanych

n j n

j

G 



^1



Twierdzenie powyższe mówi nam, że aby nasz schemat nie był rozbieżny, to wszystkie wartości własne

G

muszą leżeć na płaszczyźnie zespolonej (bo

G

mogą być zespolone) wewnątrz okręgu o promieniu

1

.

W naszym przypadku wartości

G

były rzeczywiste, więc weźmy ciekawszy przypadek:

 0



 





x u t

u

) , ( t x x f

u t

u 



 





Uwaga: stabilność każdego problemu niehomogenicznego

badamy upraszczając równanie do postaci homogenicznej.

Zastosujmy schemat:

1

0

1

 

 





_



x u u

t u

u

_j^p _j^{p p}_j ^p_j

różnica zwykła różnica wsteczna

 

j^p

p j p

j p

j

r   



^1

 

_1



Zastępując

u





0 1 1

j n n

j

G 



^



^

Pamiętając, że

oraz, że



⁰_j

 e

^ikmx

otrzymujemy



^{p ik x p ik x x}



^{p ik x}

x

p

e

ik_m

r G e

_m

G e

_m

G e

_m

G

^1

 

^{(  )}



Dzieląc obustronnie przez

G

^p

e

^ikmx

 ^e

^{ik x}

 ^{r re}

^{ik x}

^{r re}

^i

r

G  1 

^{ m}

 1  1  

^{ m}

 1  

^

Wniosek: schemat jest niestabilny dla każdego

r.

 

j^p p

j p j p

j r   

 ^1   _1 

0 1 1

j n n

j G





^  ^

x ik

j

 e

^m



0

Życiowy przykład…

n=1000;

dh=1.0/n;

dt=1.0*dh;

for(x=0;x<=n;x++) { X=(double)x/n;

u[x][0]=0.2*exp(-200*(X-0.5)*(X-0.5));

u[x][1]=0.2*exp(-200*(X-0.5-dt)*(X-0.5-dt));

};

Warunki początkowe

for(x=0;x<=n;x++) { L=x-1;

P=x+1;

if(x==0) L=n;

if(x==n) P=0;

u[x][2]=(u[P][1]-2*u[x][1]+u[L][1])*dt*dt/dh/dh+2*u[x][1]-u[x][0];

};

for(x=0;x<=n;x++) { u[x][0]=u[x][1];

u[x][1]=u[x][2];

};

Warunki brzegowe (periodyczne) Jeden krok czasowy:

𝜕²𝑢

𝜕𝑡² = 𝜕²𝑢

𝜕𝑥^{2 2}

1 1

2

1

1

2 2

x

u u

u t

u u

u

_i^t _i^t _i^t _i^t _i^t _i^t





 







^ _{ }



0.2 0.4 0.6 0.8 1.0 0.05

0.10 0.15 0.20

Można obliczyć warunek stabilności dla podanego schematu, ale nie jest prosto…

 ^{1 2 sin ( )}  ^{1 0}

2

^{2 2} ₂

2

  r G  

G

^

Czyli mamy równanie kwadratowe.

Schemat niejawny Eulera

- punkt uwzględniany przy obliczaniu różnicy czasowej

- punkt uwzględniany przy obliczaniu różnicy przestrzennej

Schemat jawny Schemat niejawny

Wróćmy do równania ciepła

2

1 1

1

2

x

u u

u t

u

u

_i^j _i^j _i^j _i^j _i^j





 





_{ }



Było:

2

1 1 1

1 1

1

2

x

u u

u t

u

u

_i^j _i^j _i^j _i^j _i^j





 





_^{ } _^



Teraz:

j i j

i j

i j

i

r u ru u

ru    



_^₁¹

( 1 2 )

^1 _^₁¹

 













 



















 













 













 

 

 





 

 

 



























 











1 1

3 2

1 0 1

1 1 1 3

1 2

1 1

2 1 2

1 2

1 2

1

j n j

n

j j

j j

j n j j j

ru u

u u

ru u

u u u u

r r

r r

r r

r r

r

r r



 



2

1 1 1

1 1

1

2

x

u u

u t

u

u

_i^j _i^j _i^j _i^j _i^j





 





_^{ } _^



Zatem mamy układ równań z macierzą trójdiagonalną.

który musimy rozwiązać w każdym kroku czasowym.

Stabilność schematu niejawnego

Postępując analogicznie do poprzednich przypadków przyjmujemy po paru przekształceniach

) (

sin 4

1

1

2 2 k_m x

G r

_

 

co daje ograniczenie na

G

:

4 1 1

1  

 G

r

co jest spełnione dla dowolnego

r

(oczywiście większego od zera).

Zatem schemat niejawny jest zawsze stabilny.

Schemat Cranka-Nicholsona

2

1 1 1

1 1 2

1 1

1

2

) 1

2 (

x

u u

u x

u u

u t

u

u

_i^j _i^j _i^j _i^j _i^j _i^j _i^j _i^j





 

 



 





_{  }^{ } _^



 

Możemy wyobrazić sobie bardziej ogólny schemat:

gdzie 0 ≤  ≤ 1.

Przyjmijmy  = ½.

Pierwsza pochodna czasowa jako różnica centralna w t^l+1/2

t u u

t

u

_i^l _i^l



 





^1

Druga pochodna przestrzenna jako różnica centralna ważona

 

 









 





 





_{  }^{ } _^

2

1 1 1

1 1

2 1

1 2

2

) (

2 )

( 2 2

1

x

u u

u x

u u

u x

u

_i^l _i^l _i^l _i^l _i^l _i^l

l i l

i l

i l

i l

i l

i

r u ru ru r u ru

ru

_^₁¹

 2 ( 1  )

^1



_^₁¹



_1

 2 ( 1  ) 

_1



Po paru przekształceniach

A zatem znowu mamy układ równań z macierzą trójdiagonalną.

Porównanie schematów:

Można pokazać, że układ jest zawsze stabilny.

) 2 / ( sin 2 1

) 2 / ( sin 2 1

2 2







 

G 2

3

2 2

1 0.5 0.5 1

Przewaga nad schematem niejawnym – błąd O(t²) zamiast O(t).

𝑥 𝑡 = 𝑥₀ + 𝜀(𝑡)

𝜕²𝑥

𝜕𝑡² = 𝐹 𝑥 = −𝜕𝑉

𝜕𝑥

𝑉 𝑥 = 𝑉 𝑥₀ + 𝜕𝑉

𝜕𝑥 𝑥₀𝜀 + 1 2

𝜕²𝑉

𝜕𝑥² 𝑥₀𝜀²

𝜕²𝜀

𝜕𝑡² = −𝐶𝜀 𝜀~𝑒^𝐶′𝑡

x₀ x V