Metody numeryczne rozwiązywania równań różniczkowych

Metody numeryczne rozwiązywania równań różniczkowych

Marcin Orchel

Spis treści

1 Wstęp 1

1.1 Metody przybliżone dla równań pierwszego rzędu . . . 1

1.1.1 Metoda kolejnych przybliżeń Picarda . . . 1

1.1.2 Całkowanie przez rozwinięcie w szereg . . . 3

1.2 Metody przybliżone dla równań drugiego rzędu . . . 4

1.3 Rozwiązywanie zagadnień brzegowych metodami numerycznymi . . . 6

1.3.1 Metoda różnicowa . . . 6

1.3.2 Metoda postulowania postaci rozwiązania . . . 7

2 Zadania 9 2.1 Zadania na 3.0 . . . 9

2.2 Zadania na 4.0 . . . 11

2.3 Zadania na 5.0 . . . 11

1 Wstęp

1.1 Metody przybliżone dla równań pierwszego rzędu 1.1.1 Metoda kolejnych przybliżeń Picarda

Całkując równanie różniczkowe

y⁰= f (x, y) (1)

z warunkiem początkowym y = y₀ dla x = x₀ dostajemy wzór:

y = y0+ Z x

x0

f (x, y) dx (2)

Wzór iteracyjny:

yi= y₀+ Z x

x0

f (x, yi−1) dx (3)

Przykład 1.

y⁰ = e^x− y² (4)

z warunkiem początkowym: x₀ = 0, y₀ = 0.

y₀ = 0 (5)

y1= 0 + Z x

0

e^x− y₀²dx = e^x− 1 (6)

y2 = 0 + Z x

0

e^x− (e^x− 1)²dx = Z x

0

3e^x− e^2x− 1dx = 3e^x− 3 −1

2e^2x+ 1

2− x (7) ...

Przykład 2.

y⁰ = 1 4x −1

4y + 2 (8)

z warunkiem początkowym: x₀ = 0, y₀ = 0:

y₀ = 0 (9)

y₁= 0 + Z _x

0

1

4x + 2dx = 2x + 1

8x² (10)

y₂= 0 + Z _x

0

1 4x −1

4

 2x +1

8x²



+ 2dx = Z _x

0

−1 4x − 1

32x²+ 2dx (11) y₂ = 2x −1

8x²− 1

96x³ (12)

y₃= Z x

0

1 4x −1

4

 2x −1

8x²− 1 96x³



+ 2dx = Z x

0

−1 4x + 1

32x²+ 1

384x³+ 2dx (13) y₃= 2x − 1

8x²+ 1

96x³+ 1

1536x⁴ (14)

... Ponieważ

e^x= 1 + x 1!+ x²

2! + . . . (15)

e^−14x0 = −1

4e^−14x (16)

e^−14x00= − 1

16e^−14x (17)

... Szereg Maclaurina ma postać

∞

X

n=0

1

n!f⁽ⁿ⁾(0) xⁿ (18)

gdzie f⁽⁰⁾(0) = f (0). Z szeregu Maclaurina otrzymujemy, że:

e^−14x= 1 −1 4x + 1

32x²− 1

384x³ (19)

x + 4 − 4e^−14x= x + 4 − 4 + x −1

8x²+ 1

96x³− 1

1536x⁴. . . (20) Rozwiązanie na wolframalpha.com,http: // www. wolframalpha. com/ input/ ?i= y% 27+

%3D+ 1% 2F4x+ -+1% 2F4y+ %2B+ 2.

1.1.2 Całkowanie przez rozwinięcie w szereg Szereg Taylora dla funkcji:

y = y (x0) + (x − x₀) y⁰(x₀) +(x − x₀)²

2 y⁰⁰(x₀) + . . . +(x − x₀)ⁿ

n! y⁽ⁿ⁾(x₀) + . . . (21) lub inaczej zapisany

∞

X

n=0

1

n!f⁽ⁿ⁾(x₀) (x − x₀)ⁿ (22) gdzie f⁽⁰⁾(x₀) = f (x₀).

Przykład 3.

y⁰ = e^x− y² (23)

Warunek początkowy: (0, 0) 1 sposób. Postulujemy rozwiązanie:

y = a₁x + a₂x²+ a₃x³+ . . . + a_nxⁿ+ . . . (24) Kwadrat szeregu potęgowego jest równy

∞

X

n=1

anxⁿ

!2

=

∞

X

n=2 n−1

X

i=1

aian−i

!

xⁿ (25)

Po podstawieniu tego szeregu do równania wyjściowego i zastosowaniu wzoru na kwadrat szeregu potęgowego otrzymujemy:

a₁+2a₂x+3a₃x²+. . .+^a²₁x²+ 2a₁a₂x³+^a²₂+ 2a₁a₃^x⁴+ . . .^= 1+x+x² 2 +x³

6 +. . . (26) Mamy następujące równania na współczynniki:

a1 = 1 (27)

2a₂ = 1 (28)

3a₃+ a²₁ = 1

2 (29)

4a₄+ 2a₁a₂ = 1

6 (30)

..., w rezultacie otrzymujemy:

y = x +x² 2 −x³

6 − 5

24x⁴+ . . . (31)

2 sposób. Obliczamy kolejne pochodne:

y⁰⁰ = e^x− 2yy⁰ (32)

y⁰⁰⁰ = e^x− 2y⁰²− 2yy⁰⁰ (33) y⁽⁴⁾= e^x− 6y⁰y⁰⁰− 2yy⁰⁰⁰ (34) Wiemy, że wartość szukanej funkcji w punkcie x₀ jest y₀. Dlatego też

y⁰(x₀) = e^x0− y²₀ (35)

Z twierdzenia Taylora dla (0, 0) otrzymujemy:

y = x +x² 2 −x³

3! − 5

4!x⁴+ . . . (36)

1.2 Metody przybliżone dla równań drugiego rzędu

Jeśli znamy jedno rozwiązanie y₁ to drugie rozwiązanie szczególne równania liniowego jednorodnego wynosi:

y₂ = Ay₁ Z e⁻

Rpdx

y₁² dx (37)

Rozwiązanie równania typu:

s (x) y⁰⁰+ p (x) y⁰+ q (x) y = F (x) (38) gdzie funkcje s (x), p (x), q (x) i F (x) są wielomianami lub w pewnym ustalonym obsza- rze dają się rozwinąć w szereg potęgowy względem x − x₀ do nich zbieżny. Rozwiązanie wyznaczamy ze pomocą metody współczynników nieoznaczonych, rozwiązanie postulu- jemy w postaci

y = a₀+ a₁(x − x₀) + a₂(x − x₀)²+ . . . (39) i podstawiamy do równania wyjściowego. Przyrównując współczynniki przy x-ach otrzy- mujemy wartości parametrów a_i.

Przykład 4.

y⁰⁰+ xy = 0 (40)

Podstawiamy do niego

y = a0+ a₁x + a2x²+ a₃x³+ a₄x⁴+ a₅x⁵. . . (41) y⁰= a₁+ 2a₂x + 3a₃x²+ 4a₄x³+ 5a₅x⁴. . . (42) y⁰⁰= 2a₂+ 6a₃x + 12a₄x²+ 20a₅x³. . . (43) i otrzymujemy:

2a₂+ 6a₃x + 12a₄x²+ 20a₅x³+ . . . + a₀x + a₁x²+ a₂x³+ a₃x⁴+ a₄x⁵+ a₅x⁶= 0 (44)

2a₂ = 0 6a₃+ a₀ = 0 12a₄+ a₁ = 0 20a₅+ a₂ = 0 . . .

(45)

Rozwiązując je otrzymujemy:

a2 = 0 a₃ = −_2·3^a0 a₄ = −_3·4^a1 a5 = 0 . . .

(46)

Wszystkie współczynniki są 0 lub zależne od a₀ i a₁. Po ich podstawieniu do y otrzymamy rozwiązanie ogólne postaci:

y = a₀ 1 − x³

2 · 3+ x⁶

2 · 3 · 5 · 6− . . .

!

+ a₁ x − x⁴

3 · 4 + x⁷

3 · 4 · 6 · 7− . . .

!

(47) Równanie na wolframalpha.com, http: // www. wolframalpha. com/ input/ ?i= y% 27%

27% 2Bxy% 3D0.

Rozwiązanie równania typu

x²y⁰⁰+ xp (x) y⁰+ q (x) y = 0 (48) Rozwiązujemy za pomocą metody współczynników nieoznaczonych, przy założeniu, że p (x) i q (x) dają się rozwinąć w szereg potęgowy względem x. Rozwiązania wtedy mają postać:

y = x^ra0+ a₁x + a2x²+ . . .^ (49) gdzie r jest pierwiastkiem równania wskaźnikowego (określającego)

r (r − 1) + p (0) r + q (0) = 0 (50)

Jeśli pierwiastki równania wskaźnikowego są różne i ich różnica nie jest liczbą całkowi- tą, to otrzymujemy w ten sposób dwa niezależne liniowo rozwiązania. W przeciwnym razie metoda współczynników nieoznaczonych daje tylko jedno rozwiązanie. Drugie roz- wiązanie możemy wtedy wyznaczyć ze wzoru (37). Przykładem tego typu równania jest równanie Bessela postaci:

x²y⁰⁰+ xy⁰+^x²− n^2y = 0 (51) Przykład 5.

x²y⁰⁰+ xy⁰+^x²− 4^y = 0 (52) Rozwiązanie na wolframalpha.com,http: // www. wolframalpha. com/ input/ ?i= x^2y%

27% 27+ %2B+ xy% 27+ %2B+ %28x^2-4% 29y+ %3D+ 0.

1.3 Rozwiązywanie zagadnień brzegowych metodami numerycznymi 1.3.1 Metoda różnicowa

Po angielsku finite difference method. Metoda ta będzie przedstawiona na przykładzie zagadnienia brzegowego dla równania liniowego drugiego rzędu:

y⁰⁰(x) + p (x) y⁰(x) + q (x) y (x) = f (x) (53) gdzie

a ≤ x ≤ b (54)

y (a) = α (55)

y (b) = β (56)

Przedział [a, b] dzielimy za pomocą położonych w równych odstępach węzłów

x_v = x₀+ vh (57)

gdzie

v = 0, 1, . . . , n; x0 = a; x_n= b (58) Zatem mamy n + 1 węzłów. Tworzymy n − 1 równań postaci:

y⁰⁰(x_v) + p (x_v) y⁰(x_v) + q (x_v) y (x_v) = f (x_v) (59) dla

v = 1, . . . , n − 1 (60)

Następnie wartości pochodnych zastępujemy wyrażeniami skończonymi (finite differen- ce):

y⁰(x_v) ≈ yv+1− y_v−1

2h (61)

y⁰⁰(x_v) ≈ yv+1− 2y_v+ y_v−1

h² (62)

Drugie wyrażenie można wyprowadzić z zastosowania do pierwszych pochodnych forward difference a do drugich backward difference: po pierwszych pochodnych mamy

y⁰(x_v−1) ≈ y_v− y_v−1

h (63)

y⁰(x_v) ≈ y_v+1− y_v

h (64)

i wreszcie

y⁰⁰(x_v) ≈ y_v⁰ − y_v−1⁰

h (65)

Po podstawieniu otrzymujemy układ n − 1 równań liniowych, gdzie niewiadomymi są wartości funkcji w węzłach, y_v. Wartości w pierwszym i ostatnim węźle są znane.

1.3.2 Metoda postulowania postaci rozwiązania

Metoda ta będzie przedstawiona na przykładzie zagadnienia brzegowego dla równania liniowego drugiego rzędu:

y⁰⁰(x) + p (x) y⁰(x) + q (x) y (x) = f (x) (66) gdzie

a ≤ x ≤ b (67)

y (a) = α (68)

y (b) = β (69)

Postulujemy rozwiązanie postaci:

g (x) =

n

X

i=0

a_ig_i(x) (70)

gdzie każda funkcja g_i(x) spełnia warunki brzegowe i funkcje te są liniowo niezależne.

Po podstawieniu tego rozwiązania do równania wyjściowego otrzymujemy:

ε (x; a₁, a₂, . . . , a_n) = g⁰⁰(x) + p (x) g⁰(x) + q (x) g (x) − f (x) (71) gdzie e jest błędem przybliżenia, tzw. defektem. Można zauważyć, że defekt jest funkcją liniową współczynników a_i.

Często postuluje się rozwiązanie postaci:

g (x) = g₀(x) +

n

X

i=1

a_ig_i(x) (72)

gdzie g₀(x) spełnia zadane warunki brzegowe, natomiast g_i(x) warunki:

gi(a) = g_i(b) = 0 (73)

i = 1, 2, . . . , n (74)

Wtenczas dla dowolnych wartości a_i spełnione są warunki brzegowe. Przykład funkcji g₀(x):

g₀(x) = α +β − α

b − a (x − a) (75)

Widzimy, że są dla niej spełnione warunki (68), (69).

Metody wyznaczania współczynników a_i.

• Metoda uzgodnienia. Zakładamy, że błąd znika w n punktach:

ε (x_v; a₁, a₂, . . . , a_n) = 0 (76) gdzie

v = 1, 2, . . . , n (77)

a < x₁ < x₂< . . . < xn< b (78) Otrzymujemy układ n równań liniowych na współczynniki a_i.

• Metoda kwadratu błędu. Minimalizujemy:

F (a₁, a₂, . . . , a_n) = Z b

a

ε²(x; a₁, a₂, . . . , a_n) dx (79) Warunek konieczny na ekstremum dostarcza nam układ równań liniowych:

∂F

∂ai

= 0 (80)

dla

i = 1, 2, . . . , n (81)

• Metoda Galerkina. Wymuszamy ortogonalność błędu, tzn. musi zachodzić:

Z b a

 (x; a₁, a₂, . . . , a_n) g_i(x) dx = 0 (82) dla i = 1, 2, . . . , n, otrzymujemy w ten sposób układ równań liniowych na współ- czynniki a_i.

• Metoda Ritza. Często rozwiązanie zagadnienia brzegowego jest równocześnie roz- wiązaniem pewnego zagadnienia wariacyjnego, tzn. y (x) minimalizuje całkę po- staci:

I [y] = Z b

a

H x, y, y⁰
dx (83)

Jeśli znamy funkcję H (x, y, y⁰) to wstawiamy g (x) za y (x). Następnie minimali- zujemy wyrażenie

I [y] = I (a₁, a₂, . . . , a_n) (84) Z warunków koniecznych na ekstremum:

∂I

∂a_i = 0 (85)

dla

i = 1, 2, . . . , n (86)

otrzymujemy n równań na współczynniki a_i.

Przykład 6. Pod ściśle określonymi warunkami nałożonymi na funkcje p, q, f i y zagadnienie brzegowe:

−^p (x) y⁰(x)^0+ q (x) y (x) = f (x) (87)

y (a) = α (88)

y (b) = β (89)

jest równoważne zagadnieniu wariacyjnemu:

miny I (y) = Z b

a

hp (x) y⁰²(x) + q (x) y²(x) − 2f (x) y (x)ⁱdx (90)

y (a) = α (91)

y (b) = β (92)

2 Zadania

2.1 Zadania na 3.0

• Rozwiązać metodą numeryczną dowolne równanie zwyczajne z warunkiem począt- kowym oraz zwyczajne drugiego rzędu z warunkami brzegowymi oraz cząstkowe z wybranymi warunkami

• wyświetlić wykres fazowy, dla różnych wartości t wyświetlamy punkty y(t) i y⁰(t)

• rozwiązać równanie opisujące oscylator van der Pola

y⁰⁰(t) + ε^y²(t) − 1^y⁰(t) + y (t) = 0 (93) dla t ∈ (0, t_max) z warunkiem początkowym y (0) = 0, 01, y⁰(0) = 0.

Równanie to można przekształcić do układu równań pierwszego rzędu

y⁰₁ = y₂ (94)

y⁰₂= −ε^y²₁− 1^y2− y₁ (95) Funkcja y₁ odpowiada y, y₂ odpowiada y⁰.

Naszkicować wykres y(t) oraz wykres fazowy dla t ∈ [0, 60], ε = 0, 8. Zaobserwować dokładność wykresu fazowego dla różnych wartości kroku.

• rozwiązać równanie van der Pola dla ε = 0 oraz dla dużego t_max, np. t_max = 0, 2 · 10⁵. Na wykresie y(t) widzimy wygaszaną funkcję tymczasem norma jest stała i równa 0,01

• Rozwiązać układ równań różniczkowych zwyczajnych Lorentza dy1

dt = −σ (y₁− y₂) (96)

dy₂

dt = −y₁y₃+ ry₁− y₂ (97)

dy₃

dt = y₁y₂− by₃ (98)

dla σ = 10, r = 28, b = 8/3. Aby sprawdzić czy rozwiązania są stabilne, porównać rozwiązania dla różnych tolerancji błędu. Inne wartości do przetestowania to r = 99, 96. Naszkicować wykres o współrzędnych (y₁, y₂, y₃) dla różnych wartości t.

• Porównać z rozwiązaniem symbolicznym z wolframalpha.com.

Wskazówki

•

Wskazówki do R

• pakiet deSolve

• równanie van der Pola

• https://www.rdocumentation.org/packages/deSolve/topics/ode

• https://www.rdocumentation.org/packages/deSolve/topics/lsode

• https://journal.r-project.org/archive/2010/RJ-2010-013/RJ-2010-013.pdf Wskazówki do Matlaba

• Dla równań różniczkowych cząstkowych brak rozwiązania symbolicznego w Matla- bie.

• http://www.mathworks.com/help/matlab/ref/pdepe.html

• http://www.mathworks.com/help/matlab/ordinary-differential-equations.

html

• http://www.mathworks.com/help/matlab/boundary-value-problems.html

• http://www.mathworks.com/help/matlab/partial-differential-equations.

html

Przykłady w R

• library(deSolve)

## Time

t <- seq(0, 100, 1)

## Initial population N0 <- 10

## Parameter values

params <- list(r=0.1, K=1000)

## The logistic equation

fn <- function(t, N, params) with(params, list(r * N * (1 - N / K)))

## Solving and plotin the solution numerically out <- ode(N0, t, fn, params)

plot(out, lwd=2, main="Logistic equation\nr=0.1, K=1000, N0=10")

## Ploting the analytical solution

with(params, lines(t, K * N0 * exp(r * t) / (K + N0 * (exp(r * t) - 1)), col=2, lwd=2))

• równanie pierwszego rzędu library(deSolve)

t <- seq(0, 100, 1) N0 <- 10

fn <- function(t, N, parms) list(N * (1 - N)) out <- ode(N0, t, fn)

plot(out, lwd=2, main="")

• równanie drugiego rzędu zamienione na układ równań pierwszego rzędu t <- seq(0, 100, 1)

N0 <- c(10, 10)

fn <- function(t, N, parms) list(c(y[2], -y[1]-t)) out <- ode(N0, t, fn)

plot(out, lwd=2, main="")

2.2 Zadania na 4.0

Rozwiązać metodą numeryczną następujące zagadnienia brzegowe:

y⁰⁰(x) + y (x) = 0 (99)

z warunkami:

y (0) = 0 (100)

y

π 2



= 2 (101)

Rozwiązać również powyższe równania symbolicznie na wolframalpha.com.

2.3 Zadania na 5.0

Rozwiązać metodą numeryczną następujący problem graniczny:

∂u

∂t = ∂²u

∂x² (102)

z warunkami:

u (x, 0) = sin πx (103)

dla

x ∈ [0, 1] (104)

oraz

u (0, 1) = u (1, t) = 0 (105)

Rozwiązać również powyższe równania symbolicznie na wolframalpha.com.

Literatura

[1] I. N. Bronsztejn, K. Siemiendiajew, G. Musiol, and H. Möhlig, Nowoczesne kompen- dium matematyki. Wydawnictwo naukowe PWN, 2004.

[2] J. Niedoba and W. Niedoba, Równania różniczkowe zwyczajne i cząstkowe. Wydaw- nictwa AGH, 2001.