Metody numeryczne

Równania ró˙zniczkowe zwyczajne, cz˛e´s´c II

Janusz Szwabi´nski

[e-mail]

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.1/33

Równania ró˙zniczkowe zwyczajne, cz˛e´s´c II

1. Zagadnienia pocz ˛atkowe (ci ˛ag dalszy)

•

Metody niejawne

•

Algorytm Verleta

2. Zagadnienia brzegowe i własne

•

Prosta metoda strzałów

•

Równania liniowe

•

Jednowymiarowe równanie Schrödingera

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.2/33

Metody niejawne (1)

y

= (y

1

, . . . , y

n

), g(y, t) = [g

1

(y, t), . . . , g

n

(y, t)]

dy

dt

= g(y, t), y(t

0

) = y

0

Formalne rozwi ˛azanie

y

n

+m

= y

n

+

Z

t

n

+m

t

n

g

(y, t)dt

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.3/33

Metody niejawne (2)

m

= 1

i metoda trapezów:

y

n

+1

' y

n

+

τ

2

[g(y

n

, t

n

) + g(y

n

+1

, t

n

+1

)]

⇒ równanie wzgl˛edem y

n

+1

0 = y

n

− y

n

+1

+ +

τ

2

[g(y

n

, t

n

) + g(y

n

+1

, t

n

+1

)]

Metoda parabol

⇒ wzór Adamsa–Moultona

y

n

+1

= y

n

+

τ

15

[5g(y

n

+1

, t

n

+1

)+8g(y

n

, t

n

)−g(y

n

−1

, t

n

−1

)]+O(τ

4

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.4/33

Metody niejawne (3)

Własno´sci

•

du˙zy nakład oblicze´n

•

stabilno´s´c

Je´sli g(y, t) jest liniow ˛a funkcj ˛a y

g

(y, t) = ˜

g

(t)y

⇒ równanie mo˙zna rozwi ˛aza´c analitycznie

y

n

+1

=



1 +

1
2

˜

g

(t

n

)τ

1 −

1
2

˜

g

(t

n

+1

)τ



y

n

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.5/33

Algorytm Verleta (1)

¨

~r

i

=

1

m

i

~

F

i

~r

i

- wektor poło˙zenia i–tej cz ˛astki

m

i

-

masa i–tej cz ˛astki

~

F

i

- siła działaj ˛aca na i–t ˛a cz ˛astk˛e

Je´sli liczba cz ˛astek jest du˙za

⇒ układ wielu równa´n ró˙zniczkowych

⇒ wyliczenie siły mo˙ze by´c skomplikowane

⇒ algorytm Rungego–Kutty raczej si˛e nie nadaje ze wzgl˛edu na

wydajno´s´c

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.6/33

Algorytm Verleta (2)

Rozwijamy funkcje ~r

i

(t + ∆t)

i ~r

i

(t − ∆t) w szereg Taylora

wokół punktu t:

~r

i

(t + ∆t) = ~r

i

(t) + ∆t~v

i

(t) +

1
2

∆t

2

~a

i

(t) + . . .

~r

i

(t − ∆t) = ~r

i

(t) − ∆t~v

i

(t) +

1
2

∆t

2

~a

i

(t) + . . .

Dodaj ˛ac stronami, otrzymamy

~r

i

(t + ∆t) = 2~r

i

(t) − ~r

i

(t − ∆t) + ∆t

2

~

F

i

(t)

m

i

+ O(∆t

4

)

Odejmowanie stronami prowadzi do

~v

i

(t) =

~r

i

(t + ∆t) − ~r

i

(t − ∆t)

2∆t

+ O(∆t

3

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.7/33

Algorytm Verleta (3)

Zalety

•

jedno wyliczenie siły w ka˙zdym kroku czasowym

•

du˙za dokładno´s´c

•

stabilno´s´c

Wady

•

algorytm nie startuje z jednego punktu

•

pr˛edko´s´c nie jest generowana automatycznie, trzeba j ˛a
dodatkowo wylicza´c

•

pr˛edko´s´c obarczona jest zwykle wi˛ekszym bł˛edem, ni˙z
poło˙zenie

⇒ do wyliczenia pr˛edko´sci u˙zyj zmodyfikowanych wersji tego

algorytmu

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.8/33

Zagadnienia brzegowe

Typowe fizyczne zagadnienie brzegowe ma posta´c równania
ró˙zniczkowego drugiego rz˛edu

u

00

(x) = f (u, u

0

, x

)

z zadanymi na brzegach warto´sciami u lub u

0

W 1D mamy cztery rodzaje warunków brzegowych:

1. u(0) = u

0

i u(1) = u

1

2. u(0) = u

0

i u

0

(1) = v

1

3. u

0

(0) = v

0

i u(1) = u

1

4. u

0

(0) = v

0

i u

0

(1) = v

1

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.9/33

Zagadnienia własne

Je´sli

•

równanie ró˙zniczkowe zawiera parametr

•

rozwi ˛azanie równania spełnia zagadnienie brzegowe tylko
dla pewnych warto´sci parametru

•

musimy ten parametr wyznaczy´c

⇒ zagadnienie własne

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.10/33

Prosta metoda strzałów (1)

Rozwa˙zmy zagadnienie brzegowe

u

00

= f (u, u

0

, x

), u(0) = α, u(1) = β

Zagadnienie pocz ˛atkowe

u

00

= f (u, u

0

, x

), u(0) = α, u

0

(0) = s

ma wówczas na ogół dokładnie jedno rozwi ˛azanie

u

(x) ≡ u(x, s)

zale˙z ˛ace od wyboru warto´sci pocz ˛atkowej s dla u

0

(0)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.11/33

Prosta metoda strzałów (2)

Szukamy takiego s = ˜s, aby

u

(0) = u(1, ˜

s

) = β

Innymi słowy, szukamy miejsce zerowego ˜s funkcji

F

(s) ≡ u(1, s) − β

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.12/33

Prosta metoda strzałów (3)

Metoda bisekcji

u(x,s )

(1)

u(x,s )

(0)

~

u(x,s )

x

0

1

u

β

α

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.13/33

Prosta metoda strzałów (4)

Je´sli F (s) jest ró˙zniczkowalne wzgl˛edem s

⇒ metoda Newtona

s

(i+1)

= s

(i)

−

F

(s

(i)

)

F

0

(s

(i)

)

Mo˙zemy zast ˛api´c F

0

(s)

odpowiednim ilorazem ró˙znicowym

∆F (s

(i)

) =

F

(s

(i)

+ ∆s

(i)

) − F (s

(i)

)

∆s

(i)

Trudno´sci:

•

zbyt du˙ze ∆s

(i)

⇒ przybli˙zenie niezbyt dokładne

•

zbyt małe ∆s

(i)

⇒ znoszenie si˛e składników

⇒ ∆s

(i)

=

√

s

(i)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.14/33

Równania liniowe (1)

u

00

(x) + d(x)u

0

(x) + q(x)u(x) = s(x)

u

(0) = u

0

, u

(1) = u

1

Funkcje d(x), q(x) i s(x) s ˛a gładkie

⇒ zasada superpozycji rozwi ˛aza´n

⇒ nie musimy szuka´c dokładnie ˜s w metodzie strzałów

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.15/33

Równania liniowe (2)

u

(x, s

1

), u(x, s

2

)

- rozwi ˛azania próbne

Szukamy rozwi ˛azania postaci

u

(x) = au(x, s

1

) + bu(x, s

2

)

Z warunków brzegowych

a

+ b = 1

u

(1, s

1

)a + u(1, s

2

)b = u

1

czyli

a

=

u

(1, s

2

) − u

1

u

(1, s

2

) − u(1, s

1

)

, b

=

u

1

− u(1, s

1

)

u

(1, s

2

) − u(1, s

1

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.16/33

Zagadnienia Sturma–Liouville’a (1)

[p(x)u

0

(x)]

0

+ q(x)u(x) = s(x)

Przykłady

•

równanie Legendre’a

•

równanie Bessela

Ze wzorów trójpunktowych mamy:

∆

1

≡

u

n

+1

− u

n

−1

2h

= u

0
n

+

h

2

6

u

(3)
n

+ O(h

4

)

∆

2

≡

u

n

+1

+ u

n

−1

− 2u

n

h

2

= u

00
n

+

h

2

12

u

(4)
n

+ O(h

4

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.17/33

Zagadnienia Sturma–Liouville’a (2)

St ˛ad

p

0
n

∆

1

+ p

n

∆

2

= [p

n

u

0
n

]

0

+

h

2

12

(p

n

u

(4)
n

+ 2p

0
n

u

(3)
n

) + O(h

4

)

Zaniedbuj ˛ac drugi wyraz po prawej stronie otrzymamy

(2p

n

+ hp

0
n

)u

n

+1

+ (2p

n

− hp

0
n

)u

n

−1

= 4p

n

u

n

+ 2h

2

(s

n

− q

n

u

n

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.18/33

Zagadnienia Sturma–Liouville’a (3)

Przykład Równanie Legendre’a

d

dx



(1 − x

2

)

du
dx



+ l(l + 1)u = 0, x ∈ h−1, 1i

Załó˙zmy, ˙ze

•

nie znamy warto´sci l

•

znamy za to warto´sci wielomianu P

1

(x) = x

dla dwóch

pierwszych punktów

Szukamy l odpowiadaj ˛acego zadanym warunkom brzegowym

⇒ zagadnienie własne

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.19/33

Zagadnienia Sturma–Liouville’a (4)

Podprogram do całkowania równania Legendre’a najprostszym

algorytmem Sturma–Liouville

SUBROUTINE SLL (N,H,L,U)

IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I
DOUBLE PRECISION, INTENT (IN) :: H,L
DOUBLE PRECISION :: H2,Q,X,P,P1
DOUBLE PRECISION, INTENT (OUT), DIMENSION (N) :: U

!

H2 = 2.0*H*H
Q = L*(1.0+L)
DO I = 2, N-1

X

=

(I-1)*H-1.0

P

=

2.0*(1.0-X*X)

P1 = -2.0*X*H
U(I+1) = ((2.0*P-H2*Q)*U(I)+(P1-P)*U(I-1))/(P1+P)

END DO

END SUBROUTINE SLL

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.20/33

Zagadnienia Sturma–Liouville’a (5)

Zró˙zniczkujmy równanie Sturma–Liouville’a dwukrotnie

pu

(4)

+ 2p

0

u

(3)

= s

00

− 3p

00

u

00

− p

(3)

u

0

− p

0

u

(3)

− q

00

u

− 2q

0

u

0

− qu

00

Jednokrotne ró˙zniczkowanie daje

u

(3)

=

1
p

(s

0

− p

00

u

0

− 2p

0

u

00

− q

0

u

− qu

0

)

St ˛ad

c

n

+1

u

n

+1

+ c

n

−1

u

n

−1

= c

n

u

n

+ d

n

+ O(h

6

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.21/33

Zagadnienia Sturma–Liouville’a (6)

gdzie

c

n

+1

= 24p

n

+ 12hp

0
n

+ 2h

2

q

n

+ 6h

2

p

00
n

− 4h

2

(p

0

n

)

2

p

n

+ h

3

p

(3)
n

+2h

3

q

0

n

− h

3

q

n

p

0

n

p

n

− h

3

p

00
n

p

0

n

p

n

c

n

= 48p

n

− 20h

2

q

n

− 8h

2

(p

0

n

)

2

p

n

+ 12h

2

p

00
n

+ 2h

4

q

0

n

p

0

n

p

n

− 2h

4

q

00

n

c

n

−1

= 24p

n

− 12hp

0
n

+ 2h

2

q

n

+ 6h

2

p

00
n

− 4h

2

(p

0

n

)

2

p

n

− h

3

p

(3)
n

−2h

3

q

0

n

+ h

3

q

n

p

0

n

p

n

+ h

3

p

00
n

p

0

n

p

n

d

n

= 24h

2

s

n

+ 2h

4

s

00
n

− 2h

4

s

0
n

p

0

n

p

n

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.22/33

Algorytm Numerowa (1)

u

00

+ q(x)u = s(x)

Ze wzoru trójpunktowego mamy

u

n

+1

− 2u

n

+ u

n

−1

h

2

= u

00
n

+

h

2

12

u

(4)
n

+ O(h

4

)

Z drugiej strony

u

(4)
n

=

d

2

dx

2

(s(x) − q(x)u(x))|

x

=x

n

=

s

n

+1

− 2s

n

+ s

n

−1

h

2

−

q

n

+1

u

n

+1

− 2q

n

u

n

+ q

n

−1

u

n

−1

h

2

+ O(h

2

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.23/33

Algorytm Numerowa (2)

St ˛ad

c

n

+1

u

n

+1

+ c

n

−1

u

n

−1

= c

n

u

n

+ d

n

+ O(h

6

)

gdzie

c

n

+1

= 1 +

h

2

12

q

n

+1

c

n

= 2 −

5
6

h

2

q

n

c

n

−1

= 1 +

h

2

12

q

n

−1

d

n

=

h

2

12

(s

n

+1

+ 10s

n

+ s

n

−1

)

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.24/33

Algorytm Numerowa (3)

SUBROUTINE NUMEROV (N,H,Q,S,U)

IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I
DOUBLE PRECISION,INTENT (IN) :: H
DOUBLE PRECISION :: G,C0,C1,C2,D,UTMP
DOUBLE PRECISION, INTENT (IN), DIMENSION (N)

:: Q,S

DOUBLE PRECISION, INTENT (INOUT), DIMENSION (N) :: U

!

G = H*H/12.0

!

DO I = 2, N-1

C0 = 1.0+G*Q(I-1)
C1 = 2.0-10.0*G*Q(I)
C2 = 1.0+G*Q(I+1)
D

= G*(S(I+1)+S(I-1)+10.0*S(I))

UTMP

= C1*U(I)-C0*U(I-1)+D

U(I+1) = UTMP/C2

END DO

END SUBROUTINE NUMEROV

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.25/33

Jednowymiarowe równanie Schrödingera (1)

−

1
2

d

2

ψ

(x)

dx

2

+ V (x)ψ(x) = Eψ(x)

V

(x)

- pewien (znany) potencjał

E

-

szukana warto´s´c własna

ψ

(x)

-

szukana funkcja własna

Zachodzi

ψ

00

(x) + 2[E − V (x)]ψ(x) = 0

⇒ do całkowania równanie Schrödingera mo˙zemy u˙zy´c

algorytmu Numerowa

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.26/33

Jednowymiarowe równanie Schrödingera (2)

Szukamy stanów zwi ˛azanych

⇒ funkcja falowa powinna spełnia´c warunek

lim

|x|→∞

ψ

(x) = 0

Wady prostej metody strzałów:

•

całkowanie w kierunku wykładniczego spadku funkcji
falowej mo˙ze prowadzi´c do wzrostu bł˛edu oblicze´n

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.27/33

Jednowymiarowe równanie Schrödingera (3)

Modyfikacja metody strzałów

x

r

x

V(x)

u

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.28/33

Jednowymiarowe równanie Schrödingera (4)

Warunek ci ˛agło´sci dla obu rozwi ˛aza´n

ψ

l

(x

r

) = ψ

r

(x

r

)

ψ

0

l

(x

r

) = ψ

0

r

(x

r

)

czyli

ψ

0

l

(x

r

)

ψ

l

(x

r

)

=

ψ

0

r

(x

r

)

ψ

r

(x

r

)

Korzystaj ˛ac ze wzoru trójpunktowego otrzymamy

F

(E) =

[ψ

l

(x

r

+ h) − ψ

l

(x

r

− h)] − [ψ

r

(x

r

+ h) − ψ

r

(x

r

− h)]

2hψ

r

(x

r

)

= 0

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.29/33

Jednowymiarowe równanie Schrödingera (5)

Zarys post˛epowania:

(i) wybieramy obszar całkowania tak, aby wpływ wyboru

granic był zaniedbywalny

(ii) podajemy rozs ˛adn ˛a warto´s´c startow ˛a dla warto´sci własnej

E

(iii) wyznaczamy funkcj˛e ψ

l

(x)

, całkuj ˛ac równanie od lewej

strony do punktu x

r

+ h

, oraz funkcj˛e ψ

r

(x)

, całkuj ˛ac od

prawej strony do punktu x

r

− h. Wybieramy zero jako

pierwsz ˛a warto´s´c funkcji ψ

l

(x)

i ψ

r

(x)

(warunki

brzegowe), oraz pewn ˛a niewielk ˛a liczb˛e jako warto´s´c tych
funkcji dla drugiego punktu całkowania

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.30/33

Jednowymiarowe równanie Schrödingera (6)

(iv) aby dopasowa´c do siebie rozwi ˛azania, jedno nale˙zy

przeskalowa´c, np.:

ψ

l

(x) → ψ

l

(x)

ψ

r

(x

r

)

ψ

l

(x

r

)

(v) wyznaczamy F (E)

(vi) stosujemy jedn ˛a z metod poszukiwania pierwiastka, aby

znale´z´c warto´s´c własn ˛a z równania F (E) = 0 z wymagan ˛a
dokładno´sci ˛a

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.31/33

Jednowymiarowe równanie Schrödingera (7)

Przykład Warto´sci i funkcje własne oscylatora harmonicznego
(m = ~ = ω = 1)

E

0

E

1

E

2

E

3

E

4

0,5000000 1,5000000 2,5000000 3,5000000 4,4999999

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.32/33

Jednowymiarowe równanie Schrödingera (8)

-10

-5

0

5

10

x

-0,5

0

0,5

ψ

ψ

0

ψ

1

ψ

2

nm_slides-13.tex – Metody numeryczne – Janusz Szwabi´nski – 13/1/2003 – 22:45 – p.33/33