Почему мой код, использующий 4-й Рунге-Кутта, не дает ожидаемых значений?

У меня небольшие проблемы с попыткой понять, что не так с моим кодом, любая помощь будет чрезвычайно полезна. Я хотел решить это простое уравнение

введите здесь описание изображения

Однако значения, которые дает мой код, не совпадают с моими книжными или вольфрамовыми, поскольку y растет по мере роста x.

import matplotlib.pyplot as plt
from numpy import exp
from scipy.integrate import ode

# initial values 
y0, t0 = [1.0], 0.0

def f(t, y):
    f = [3.0*y[0] - 4.0/exp(t)]
    return f

# initialize the 4th order Runge-Kutta solver
r = ode(f).set_integrator('dopri5')
r.set_initial_value(y0, t0)

t1 = 10
dt = 0.1
x, y = [], []
while r.successful() and r.t < t1:
    x.append(r.t+dt); y.append(r.integrate(r.t+dt))
    print(r.t+dt, r.integrate(r.t+dt))

python numpy matplotlib runge-kutta
person Orbitz    schedule 18.04.2017    source источник
comment
Попробуйте уменьшить dt. Вы заметите, что ушли дальше, прежде чем он взорвется. Когда вам нужно делать слишком маленькие шаги, когда точное решение довольно гладкое, вы должны подозревать, что ваше ОДУ жесткое: en.wikipedia.org/wiki/Stiff_equation   -  person Klimaat    schedule 19.04.2017
comment
О, я вижу! С меньшим dt уравнение работает нормально, теперь я понял, я думал, что делаю что-то не так, потому что я впервые использую эту функцию оды. Ну, что ж, спасибо.   -  person Orbitz    schedule 19.04.2017
comment
@климат. Не могли бы вы опубликовать свой комментарий в качестве ответа, чтобы OP мог его принять, и вопрос можно было пометить как отвеченный? Кроме того, вы оба получите немного репутации, что приятно.   -  person Mad Physicist    schedule 19.04.2017

Ответы (1)


arrow_upward
1
arrow_downward

Ваше уравнение вообще имеет решение

y(x) = (y0-1)*exp(3*x) + exp(-x)

Из-за выбора начальных условий точное решение не содержит растущей составляющей первого слагаемого. Однако небольшие возмущения из-за дискретизации и ошибок с плавающей запятой приведут к ненулевому коэффициенту в растущем члене. Теперь в конце интервала интегрирования этот случайный коэффициент умножается на exp(3*10)=1.107e+13, что увеличит небольшие ошибки дискретизации размера 1e-7 до вклада в результат размера 1e+6, как это наблюдалось при выполнении исходного кода.

Вы можете заставить интегратор быть более точным в своих внутренних шагах, не уменьшая размер шага вывода dt, установив пороги ошибок, как в

r = ode(f).set_integrator('dopri5', atol=1e-16, rtol=1e-20)

Однако вы не можете полностью избежать ухудшения результата, так как ошибки с плавающей запятой размера 1e-16 увеличиваются до глобальных вкладов ошибок размера 1e-3.

Кроме того, вы должны заметить, что каждый вызов r.integrate(r.t+dt) будет продвигать интегратор на dt, так что сохраненный массив и напечатанные значения будут синхронизированы. Если вы хотите просто распечатать текущее состояние интегратора, используйте

    print(r.t,r.y,yexact(r.t,y0))

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

def yexact(x,y0): 
    return [ (y0[0]-1)*exp(3*x)+exp(-x) ]
person Lutz Lehmann    schedule 19.04.2017