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

(1)

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

Marcin Orchel

Spis treści

1 Wstęp 1

2 Zadania 1

2.1 Zadania na 3.0 . . . 1

1 Wstęp 2 Zadania

2.1 Zadania na 3.0

• Rozwiązać metodą numeryczną wybrane równanie różniczkowe cząstkowe pierwszego rzędu z wybranymi warunkami.

• Rozwiązać metodą numeryczną wybrane równanie różniczkowe cząstkowe drugiego rzędu z wybranymi warunkami.

• Rozwiązać metodą numeryczną wybrany układ równań różniczkowych cząstkowych pierwszego, drugiego rzędu z wybranymi warunkami.

• Rozwiązać równanie eliptyczne, paraboliczne i hiperboliczne

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

• Rozwiązać równanie reakcji-dyfuzji określone na prostokącie Ω = (0, a) × (0, b) z niewiadomą u = u(x, y, t)

∂u

∂t = β∆u + f (u) (1)

u (x, y, t) = 0 (2)

u (x, y, 0) = u0(x, y) (3)

dla (x, y) ∈ Ω

• czy zwiększenie liczby węzłów siatki polepsza rozwiązanie 1

(2)

Funkcja u₀ to wartość rozwiązania w chwili początkowej. Symbol ∆ to dwuwymiarowy operator dyfuzji, tzw. operator Laplace’a

∆u = ∂²u

∂x² +∂²u

∂y² (4)

Przykładowo f (u) = u − u³. Wskazówki

• Jak sprowadzić równanie różniczkowe cząstkowe do układu równań różniczkowych zwyczajnych? Mamy równomierną siatkę węzłów p_ij = (x_i, y_j) ∈ Ω. Przyjmiemy xi = ih_x, y_j = jh_y, gdzie h_x = a/(N_x+ 1) jest krokiem siatki w kierunku x i ana- logicznie h_y = b/(N_y+ 1). Oznaczamy u_ij(t) = u(x_i, y_j, t). Zastępujemy pochodne przestrzenne ilorazami różnicowymi

∆u_ij ≈ ∆_huij = 1

h²_x (2u_i,j− u_i+1,j − u_i−1,j) + 1

h²_y (2u_i,j− u_i,j+1− u_i,j−1) (5) przybliżamy układem równań różniczkowych zwyczajnych

d

dtU (t) = β∆_h(U (t)) + f (U (t)) (6) Tablica U jest rozmiaru N_x× N_y o elementach będących funkcjami czasu U (t) = (u_ij(t)) dla 1 ≤ i ≤ N_x, 1 ≤ j ≤ N_y, a operator ∆_hjest dyskretyzacją dwuwymiaro- wego laplasjanu. We węzłach brzegowych u jest równa zeru jak w warunkach, więc nie trzeba jej wyznaczać. Trzeba rozwiązać układ N_x· N_y równań różniczkowych zwyczajnych.

Wskazówki do R

• pakiet deSolve

• https://www.rdocumentation.org/packages/deSolve/versions/1.20/topics/

ode.1D

• https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf Przykłady w R

• Aphid <- function(t, APHIDS, parameters) { deltax <- c (0.5, rep(1, numboxes - 1), 0.5) Flux <- -D * diff(c(0, APHIDS, 0)) / deltax dAPHIDS <- -diff(Flux) / delx + APHIDS * r

# the return value list(dAPHIDS ) }

2

(3)

D <- 0.3 # m2/day diffusion rate r <- 0.01 # /day net growth rate delx <- 1 # m thickness of boxes numboxes <- 60

# distance of boxes on plant, m, 1 m intervals

Distance <- seq(from = 0.5, by = delx, length.out = numboxes)

# Initial conditions: # ind/m2 APHIDS <- rep(0, times = numboxes) APHIDS[30:31] <- 1

state <- c(APHIDS = APHIDS) # initialise state variables times <-seq(0, 200, by = 1)

print(system.time(out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")))

• drugi przykład

Aphid <- function(t, APHIDS, parameters) { deltax <- c (0.5, rep(1, numboxes - 2), 0.5) dAPHIDS <- -diff(c(0, APHIDS)) / deltax

list(dAPHIDS) }

r <- 0.01 delx <- 1 numboxes <- 60

Distance <- seq(from = 0.5, by = delx, length.out = numboxes) APHIDS <- rep(0, times = numboxes)

APHIDS[30:31] <- 1

state <- c(APHIDS = APHIDS) times <-seq(0, 200, by = 1)

image(out, method = "filled.contour", grid = Distance, xlab =

"time, days", ylab = "Distance on plant, m", main = "Aphid density on a row of plants")

• trzeci przykład

3

(4)

Aphid <- function(t, APHIDS, parameters) { deltax <- c (0.5, rep(1, numboxes-2), 0.5) dAPHIDS <- -diff(c(APHIDS, 0)) / deltax list(dAPHIDS)

}

delx <- 1 numboxes <- 50

Distance <- seq(from = 0.5, by = delx, length.out = numboxes) APHIDS <- rep(1, times = numboxes)

state <- c(APHIDS = APHIDS) state <- seq(1,50)

times <-seq(0, 10, by = 0.1)

head(out[,1:2]) plot(times, out[,1])

4