Может ли Эйлер быть лучше, чем Рунге-Кутта для некоторых функций? - PullRequest
Новая
       92

Может ли Эйлер быть лучше, чем Рунге-Кутта для некоторых функций?

0 голосов
/ 31 октября 2018

Я пытаюсь решить упражнения из «Нелинейной динамики и хаоса» Стивена Строгаца. В упражнении 2.8.3, 2.8.4 и 2.8.5 ожидается реализация метода Эйлера, улучшенного метода Эйлера и метода Рунге-Кутты (4-го порядка) соответственно для задачи начальных значений dx / dt = -x; x (0) = 1, чтобы найти x (1).

Аналитически, ответ 1 / е. И я обнаружил ошибку, полученную в каждом методе. К моему большому удивлению, в Euler я получал меньше ошибок, чем в Improved Euler и Runge-Kutta!

Мой код выглядит следующим образом. Извините за потертость. 

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

to = 0
xo = 1
tf = 1

deltaT = np.zeros([5])
errorE = np.zeros([5])
errorIE = np.zeros([5])
errorRK = np.zeros([5])


for j in range(0,5):
  n = pow(10,j)
  deltat = (tf - to)/(n)

  print ("delta t is",deltat)

  deltaT[j] = deltat

  t = np.linspace(to,tf,n)
  xE = np.zeros([n])
  xIE = np.zeros([n])
  xRK = np.zeros([n])

  xE[0] = xo
  xIE[0] = xo
  xRK[0] = xo

  for i in range (1,n):
    #Regular Euler
    xE[i] = deltat*(-xE[i-1]) + xE[i-1]

    #Improved Euler
    IEintermediate = deltat*(-xIE[i-1]) + xIE[i-1]
    xIE[i] = xIE[i-1] - deltat*(xIE[i-1] + IEintermediate)/2 

    #Runge-Kutta fourth order
    k1 = -deltat*xRK[i-1]
    k2 = -deltat*(xRK[i-1] + k1/2)
    k3 = -deltat*(xRK[i-1] + k2/2)
    k4 = -deltat*(xRK[i-1] + k3)

    xRK[i] = xRK[i-1] + (k1 + 2*k2 + 2*k3 + k4)/6

    print (deltat,xE[i],xIE[i],xRK[i])

  errorE[j] = np.exp(-1) - xE[n-1]
  errorIE[j] = np.exp(-1) - xIE[n-1]
  errorRK[j] = np.exp(-1) - xRK[n-1]

Ошибки:

Для delT = 1,0

  • Ошибка Эйлера -0,6321205588285577
  • I. Ошибка Эйлера -0,6321205588285577
  • Ошибка РК: -0,6321205588285577

Для delT = 0,1

  • Ошибка Эйлера -0.019541047828557645
  • I. Ошибка Эйлера -0.039348166379443716
  • Ошибка РК -0,03869055002863331

Для delT = 0,01

  • Эйлер -0.0018501964782845493
  • I.Euler -0.003703427083890265
  • РК -0,0036972498815148747

Для delT = 0,001

  • Эйлер -0.0001840470877806366
  • I.Euler -0.00036812480143849635
  • RK-+0,00036806344222467535

Для delT = 0,0001

  • Эйлер -1.839504510836587e-05
  • I.Euler -3.67903967520844e-05
  • РК -3.678978357835039e-05

Это законно? Если нет, то почему это происходит?

1 Ответ

0 голосов
/ 31 октября 2018

Вы выполняете только n-1 шагов интеграции с размером шага h=1/n, таким образом, вы вычисляете 

exp(-(n-1)/n)=1/e*exp(1/n)

, который имеет приблизительное значение 

1/e + 1/e*1/n

Сообщенные значения ошибок в точности равны -h/e, который является первым порядком и, таким образом, заметно искажается методом Эйлера порядка 1. Значение Эйлера, точнее 

(1-1/n)^(n-1) = exp((n-1)*(-1/n-1/(2n^2)+O(1/n^3))
              = 1/e*exp(1/(2n)+..)
              = 1/e + h/(2e) + ...

Если вы адаптируете код, чтобы сделать дополнительный шаг для достижения времени 1, вы получите правильную картину ошибки. 

delta t is  1.0
Euler          0.0             0.367879441171
imp. Euler     0.5            -0.132120558829
Runge-Kutta 4  0.375          -0.00712055882856

delta t is  0.1
Euler          0.3486784401    0.0192010010714
imp. Euler     0.368540984834 -0.00066154366211
Runge-Kutta 4  0.367879774412 -3.33241056083e-07

delta t is  0.01
Euler          0.366032341273  0.00184709989821
imp. Euler     0.367885618716 -6.17754474969e-06
Runge-Kutta 4  0.367879441202 -3.09130498977e-11

delta t is  0.001
Euler          0.367695424771  0.000184016400479
imp. Euler     0.367879502531 -6.13592486265e-08
Runge-Kutta 4  0.367879441171 -4.05231403988e-15

delta t is  0.0001
Euler          0.367861046433  1.83947385133e-05
imp. Euler     0.367879441785 -6.13176398545e-10
Runge-Kutta 4  0.367879441171 -2.6645352591e-15
...