Рунге-Кутта. Решение проблемы начального значения, которую нелегко отделить

Мы должны написать программу для численного решения следующей задачи начального значения с использованием метода Рунге-Кутты 4-го порядка. Этот алгоритм не проблема, и я могу опубликовать свое решение, когда закончу.

Проблема в том, чтобы аккуратно разделить это на что-то, что я могу поместить в Рунге-Кутта.

e^(-x') = x' −x + exp(−t^3)
x(t=0) = 1

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

Обновление: если кого-то интересует решение, код ниже. Я подумал, что это интересная проблема.

import numpy as np
import matplotlib.pyplot as plt

def Newton(fn, dfn, xp_guess, x, t, tolerance):
    iterations = 0
    value = 100.
    max_iter = 100
    xp = xp_guess
    value = fn(t, x, xp)
    while (abs(value) > tolerance and iterations < max_iter):
        xp = xp - (value / dfn(t,x,xp))
        value = fn(t,x,xp)
        iterations += 1
    root = xp
    return root

tolerance = 0.00001
x_init = 1.
tmin = 0.0
tmax = 4.0
t = tmin
n = 1
y = 0.0
xp_init = 0.5

def fn(t,x,xp):
    '''
    0 = x' - x + e^(-t^3) - e^(-x')
    '''
    return (xp - x + np.e**(-t**3.) - np.e**(-xp))

def dfn(t,x,xp):
    return 1 + np.e**(-xp)

i = 0
h = 0.0001
tarr = np.arange(tmin, tmax, h)
y = np.zeros((len(tarr)))
x = x_init
xp = xp_init
for t in tarr:
    # RK4 with values coming from Newton's method
    y[i] = x
    f1 = Newton(fn, dfn, xp, x, t, tolerance)
    K1 = h * f1
    f2 = Newton(fn, dfn, f1, x+0.5*K1, t+0.5*h, tolerance)
    K2 = h * f2
    f3 = Newton(fn, dfn, f2, x+0.5*K2, t+0.5*h, tolerance)
    K3 = h * f3
    f4 = Newton(fn, dfn, f3, x+K3, t+h, tolerance)
    K4 = h * f4
    x = x + (K1+2.*K2+2.*K3+K4)/6.
    xp = f4
    i += 1

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(tarr, y)
plt.show()

person David Rinck    schedule 02.04.2011    source источник
comment
Это производная от x (t) как в показателе экспоненты слева и в первом линейном члене справа, верно? Это немного необычно, и я немного удивлен, увидев это в домашнем задании.   -  person David Z    schedule 02.04.2011
comment
Что это. Мы все очень хорошо учились в классе, поэтому он решает некоторые сложные задачи.   -  person David Rinck    schedule 02.04.2011
comment
Возможно, вам лучше спросить о math.SE, есть ли способ преобразовать это в какое-то линейное уравнение вида y '= f (t, y). Хотя я уверен, что есть методы для применения РК к уравнению, подобному этому, я не нашел никаких доказательств их использования в Числовых рецептах, которые, как мне кажется, квалифицируют эту тему как относительно неясную ;-)   -  person David Z    schedule 02.04.2011
comment
@ david-zaslavsky: См. мой ответ о том, как применить РК к подобному уравнению.   -  person btilly    schedule 02.04.2011


Ответы (3)


Для Рунге-Кутты вам нужно только численное решение, а не аналитическое.

То есть вам нужно иметь возможность написать фрагмент кода, который принимает (x, t) и возвращает y, так что exp(-y) == y - x + exp(-t**3) с точностью до ошибки округления. Этот код может выполнять своего рода алгоритм итеративного приближения, и Рунге-Кутта будет совершенно счастлив.

Это поможет?

person btilly    schedule 02.04.2011
comment
Спасибо за вашу помощь. Я все еще не уверен. Я попробую и дам вам знать. Я думаю, мне может понадобиться использовать функцию Ламберта из-за производной в e ^ (- x '). Сдаваем это в пятницу, если вы правы, я обязательно отмечу это как таковое. - person David Rinck; 03.04.2011
comment
@svaha: Готов поспорить на серьезные деньги, что я прав, и гарантирую, что вам не нужно использовать функцию Ламберта. Но это твое домашнее задание, поэтому я не буду вдаваться в подробности. - person btilly; 03.04.2011
comment
Спасибо! Я должен был доказать это самому себе, прежде чем я смог принять ответ. - person David Rinck; 09.04.2011

Wolfram Alpha сообщает, что решение будет выглядеть так: this.

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

Также полезно знать, что такой ресурс, как Wolfram Alpha, доступен вам в любое время.

PS - Что значит быть студентом или профессором во времена Интернета, Wolfram Alpha, Google, Wikipedia и т. Д.?

person duffymo    schedule 02.04.2011
comment
Я считаю, что перед тем, как я начну, лучше понять, какой будет ответ. Совершенно верно. +1 - person fabrizioM; 03.04.2011
comment
Спасибо! Что ж, в студенческие годы это был отличный класс. Профессор - физик из Брукхейвена. Он не имеет опыта в компьютерной науке, но имеет большой опыт работы в области вычислительной биологии и моделирования Монте-Карло. Даже с учетом ресурсов, которые предоставляет Интернет, я нахожу его занятие очень сложным и полезным. - person David Rinck; 03.04.2011

Записывая K вместо x - exp (-t ^ 3), мы хотим решить exp (-y) = y - K; Я получаю y = K + W (exp (-K)), где W - функция W Ламберта, например здесь

person dmuir    schedule 05.04.2011