• Nie Znaleziono Wyników

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

N/A
N/A
Protected

Academic year: 2021

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

Copied!
12
0
0

Pełen tekst

(1)

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

y0= f (x, y) (1)

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

y = y0+ Z x

x0

f (x, y) dx (2)

Wzór iteracyjny:

yi= y0+ Z x

x0

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

Przykład 1.

y0 = ex− y2 (4)

(2)

z warunkiem początkowym: x0 = 0, y0 = 0.

y0 = 0 (5)

y1= 0 + Z x

0

ex− y02dx = ex− 1 (6)

y2 = 0 + Z x

0

ex− (ex− 1)2dx = Z x

0

3ex− e2x− 1dx = 3ex− 3 −1

2e2x+ 1

2− x (7) ...

Przykład 2.

y0 = 1 4x −1

4y + 2 (8)

z warunkiem początkowym: x0 = 0, y0 = 0:

y0 = 0 (9)

y1= 0 + Z x

0

1

4x + 2dx = 2x + 1

8x2 (10)

y2= 0 + Z x

0

1 4x −1

4

 2x +1

8x2



+ 2dx = Z x

0

−1 4x − 1

32x2+ 2dx (11) y2 = 2x −1

8x2− 1

96x3 (12)

y3= Z x

0

1 4x −1

4

 2x −1

8x2− 1 96x3



+ 2dx = Z x

0

−1 4x + 1

32x2+ 1

384x3+ 2dx (13) y3= 2x − 1

8x2+ 1

96x3+ 1

1536x4 (14)

... Ponieważ

ex= 1 + x 1!+ x2

2! + . . . (15)

e14x0 = −1

4e14x (16)

e14x00= − 1

16e14x (17)

... Szereg Maclaurina ma postać

X

n=0

1

n!f(n)(0) xn (18)

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

e14x= 1 −1 4x + 1

32x2− 1

384x3 (19)

x + 4 − 4e14x= x + 4 − 4 + x −1

8x2+ 1

96x3− 1

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

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

(3)

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

y = y (x0) + (x − x0) y0(x0) +(x − x0)2

2 y00(x0) + . . . +(x − x0)n

n! y(n)(x0) + . . . (21) lub inaczej zapisany

X

n=0

1

n!f(n)(x0) (x − x0)n (22) gdzie f(0)(x0) = f (x0).

Przykład 3.

y0 = ex− y2 (23)

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

y = a1x + a2x2+ a3x3+ . . . + anxn+ . . . (24) Kwadrat szeregu potęgowego jest równy

X

n=1

anxn

!2

=

X

n=2 n−1

X

i=1

aian−i

!

xn (25)

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

a1+2a2x+3a3x2+. . .+a21x2+ 2a1a2x3+a22+ 2a1a3x4+ . . .= 1+x+x2 2 +x3

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

a1 = 1 (27)

2a2 = 1 (28)

3a3+ a21 = 1

2 (29)

4a4+ 2a1a2 = 1

6 (30)

..., w rezultacie otrzymujemy:

y = x +x2 2 −x3

6 − 5

24x4+ . . . (31)

2 sposób. Obliczamy kolejne pochodne:

y00 = ex− 2yy0 (32)

(4)

y000 = ex− 2y02− 2yy00 (33) y(4)= ex− 6y0y00− 2yy000 (34) Wiemy, że wartość szukanej funkcji w punkcie x0 jest y0. Dlatego też

y0(x0) = ex0− y20 (35)

Z twierdzenia Taylora dla (0, 0) otrzymujemy:

y = x +x2 2 −x3

3! − 5

4!x4+ . . . (36)

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

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

y2 = Ay1 Z e

Rpdx

y12 dx (37)

Rozwiązanie równania typu:

s (x) y00+ p (x) y0+ 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 − x0 do nich zbieżny. Rozwiązanie wyznaczamy ze pomocą metody współczynników nieoznaczonych, rozwiązanie postulu- jemy w postaci

y = a0+ a1(x − x0) + a2(x − x0)2+ . . . (39) i podstawiamy do równania wyjściowego. Przyrównując współczynniki przy x-ach otrzy- mujemy wartości parametrów ai.

Przykład 4.

y00+ xy = 0 (40)

Podstawiamy do niego

y = a0+ a1x + a2x2+ a3x3+ a4x4+ a5x5. . . (41) y0= a1+ 2a2x + 3a3x2+ 4a4x3+ 5a5x4. . . (42) y00= 2a2+ 6a3x + 12a4x2+ 20a5x3. . . (43) i otrzymujemy:

2a2+ 6a3x + 12a4x2+ 20a5x3+ . . . + a0x + a1x2+ a2x3+ a3x4+ a4x5+ a5x6= 0 (44)

(5)

2a2 = 0 6a3+ a0 = 0 12a4+ a1 = 0 20a5+ a2 = 0 . . .

(45)

Rozwiązując je otrzymujemy:

a2 = 0 a3 = −2·3a0 a4 = −3·4a1 a5 = 0 . . .

(46)

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

y = a0 1 − x3

2 · 3+ x6

2 · 3 · 5 · 6− . . .

!

+ a1 x − x4

3 · 4 + x7

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

x2y00+ xp (x) y0+ 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 = xra0+ a1x + a2x2+ . . . (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:

x2y00+ xy0+x2− n2y = 0 (51) Przykład 5.

x2y00+ xy0+x2− 4y = 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.

(6)

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:

y00(x) + p (x) y0(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

xv = x0+ vh (57)

gdzie

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

y00(xv) + p (xv) y0(xv) + q (xv) y (xv) = f (xv) (59) dla

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

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

y0(xv) ≈ yv+1− yv−1

2h (61)

y00(xv) ≈ yv+1− 2yv+ yv−1

h2 (62)

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

y0(xv−1) ≈ yv− yv−1

h (63)

y0(xv) ≈ yv+1− yv

h (64)

i wreszcie

y00(xv) ≈ yv0 − yv−10

h (65)

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

(7)

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:

y00(x) + p (x) y0(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

aigi(x) (70)

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

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

ε (x; a1, a2, . . . , an) = g00(x) + p (x) g0(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 ai.

Często postuluje się rozwiązanie postaci:

g (x) = g0(x) +

n

X

i=1

aigi(x) (72)

gdzie g0(x) spełnia zadane warunki brzegowe, natomiast gi(x) warunki:

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

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

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

g0(x) = α +β − α

b − a (x − a) (75)

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

Metody wyznaczania współczynników ai.

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

ε (xv; a1, a2, . . . , an) = 0 (76) gdzie

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

a < x1 < x2< . . . < xn< b (78) Otrzymujemy układ n równań liniowych na współczynniki ai.

(8)

• Metoda kwadratu błędu. Minimalizujemy:

F (a1, a2, . . . , an) = Z b

a

ε2(x; a1, a2, . . . , an) 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; a1, a2, . . . , an) gi(x) dx = 0 (82) dla i = 1, 2, . . . , n, otrzymujemy w ten sposób układ równań liniowych na współ- czynniki ai.

• 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, y0dx (83)

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

I [y] = I (a1, a2, . . . , an) (84) Z warunków koniecznych na ekstremum:

∂I

∂ai = 0 (85)

dla

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

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

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

p (x) y0(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) y02(x) + q (x) y2(x) − 2f (x) y (x)idx (90)

y (a) = α (91)

y (b) = β (92)

(9)

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 y0(t)

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

y00(t) + εy2(t) − 1y0(t) + y (t) = 0 (93) dla t ∈ (0, tmax) z warunkiem początkowym y (0) = 0, 01, y0(0) = 0.

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

y01 = y2 (94)

y02= −εy21− 1y2− y1 (95) Funkcja y1 odpowiada y, y2 odpowiada y0.

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 tmax, np. tmax = 0, 2 · 105. 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 = −σ (y1− y2) (96)

dy2

dt = −y1y3+ ry1− y2 (97)

dy3

dt = y1y2− by3 (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 (y1, y2, y3) dla różnych wartości t.

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

Wskazówki

(10)

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")

(11)

## 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:

y00(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 = 2u

∂x2 (102)

z warunkami:

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

dla

x ∈ [0, 1] (104)

(12)

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.

Cytaty

Powiązane dokumenty

We węzłach brzegowych u jest równa zeru jak w warunkach, więc nie trzeba

W niniejszej pracy przedstawione są heurystyki dla zagadnienia przepływowego z dowolnymi maszynami równoległymi uwzględniające możliwość dzielenia partii

Zagadnienie statyczne sprzężonej termodyfuzji lepkosprężystej (por.[2j) opisane jest układem pięciu równań różniczkowo-całkowych typu splotowego, określających

określające rów nież charakter w zajem nego oddziaływ ania pola cieplnego i dyfuzyjnego oraz pola przem ieszczeń.. W ynikają one z różnych ujęć term odynam

Jak widać, wyniki otrzymane metodą Eulera i metodą Rungego-Kutty są do siebie podobne, aby jednak podobieństwo to stało się wyraźniejsze (i aby dokończyć rozwiązywania

Nazwa metody wprowadzona zostaªa przez analogi¦ do podobnej metody w dziedzinie rozwi¡zywania równa« ró»niczkowych [4].. Dla ilustracji tej metody znaleziona zostanie SORN

Autorzy wskazują, że pomimo wprowadzenia szeregu inicjatyw legislacyjnych i politycznych, mających na celu wyeliminowanie dyskryminacji oraz ułatwie- nie utrzymania

[r]