Решение задачи вариационного исчисления, множество Мандельброта

Загружаем необходимые пакеты

require(ggplot2)

## Loading required package: ggplot2

require(deSolve)

## Loading required package: deSolve

require(bvpSolve)

## Loading required package: bvpSolve Loading required package: rootSolve
## 
## Attaching package: 'bvpSolve'
## 
## The following object is masked from 'package:stats':
## 
## approx

Решаем дифференциальное уравнение

\[ y''+\cos(t)y'+ty=5 \] С начальными условиями \( y(0)=1 \), \( y'(0)=-1 \).

Сначала заменяем его на систему уравнений первого порядка: \[ y_1'=y_2 \\ y_2=5-\cos(t)y_2-ty_1 \] С начальными условиями \( y_1(0)=1 \), \( y_2(0)=-1 \).

Задаем начальные условия, само уравнение и точки, где хотим посчитать функцию:

y.init <- c(y1 = 1, y2 = -1)
eq <- function(t, y, param) {
    return(list(c(y[2], 5 - cos(t) * y[2] - t * y[1])))
}
t <- seq(from = 0, to = 10, by = 0.1)

Решаем и строим график решения:

sol <- ode(y = y.init, f = eq, times = t)
sol <- data.frame(sol)
qplot(time, y1, data = sol, geom = "point")

В R можно работать с комплексными числами:

d <- complex(re = -0.5, im = 0.2)
d

## [1] -0.5+0.2i

d * d

## [1] 0.21-0.2i

abs(d)

## [1] 0.5385

z <- 0
z <- z^2 + d

Посмотрим, что происходит при разных числах \( d \):

d.re <- seq(from = -1.8, to = 0.6, by = 0.01)
d.im <- seq(from = -1.2, to = 1.2, by = 0.01)

df <- expand.grid(dre = d.re, dim = d.im)
df$d <- complex(re = df$dre, im = df$dim)
df$z <- 0

for (k in 1:35) df$z <- df$z^2 + df$d

df$h <- exp(-abs(df$z))

qplot(dre, dim, fill = h, data = df, geom = "tile")

На вимео есть 10-ти минутное видео с очень крупным увеличением множества Мандельброта

Решим задачу вариационного исчисления

Задаем краевые условия, само уравнение и точки, где хотим посчитать функцию:

y.st <- c(y1 = -1, y2 = NA)
y.fi <- c(y1 = 3 * (exp(6) - exp(-6)), y2 = NA)

eq <- function(t, y, param) {
    return(list(c(y[2], 9 * y[1] + 4.5 * exp(3 * t))))
}
t <- seq(from = 0, to = 2, by = 0.01)

Решаем и строим график решения:

sol <- bvptwp(yini = y.st, yend = y.fi, x = t, fun = eq)
sol <- data.frame(sol)
qplot(x, y1, data = sol)