Как записать начальные условия в функции scipy.integrate.ode (или другой)?

Question

Я пытаюсь решить дифференциальное уравнение с помощью функции python scipy.integrate.odeint и сравнить его с решением mathcad.

Итак, мое уравнение равно u'' + 0.106u'+ 0.006u = 0,. Проблема, в которой я застрял, - это начальные условия, u(0)=0 and u'(1)=1. Я не понимаю, как установить u'(1)=1.

Код Python:

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

def eq(u,t):
    return [u[1], -0.106*u[1]-0.006*u[0]] #return u' and u''

time = np.linspace(0, 10) 
u0 = [0,1] # initial conditions
Z = odeint(eq,u0,time) </code>


plt.plot(time, Z)
plt.xticks(range(0,10))
plt.grid(True)
plt.xlabel('time')
plt.ylabel('u(t)')
plt.show()

Mathcad код:

u''(t) + 0.106*u'(t) +0.006*u(t) = 0
u(0) = 0
u'(1) = 1
u := Odesolve(t,10)

Диаграмма Mathcad выглядит следующим образом:

https://pp.userapi.com/c850032/v850032634/108079/He1JsQonhpk.jpg
это эталон.

И мой вывод на python:

https://pp.userapi.com/c850032/v850032634/10809c/KB_HDekc8Fk.jpg
который выглядит похоже, но ясно, что u (t) неверно.

Answer 1

Это краевая задача, вам нужно использовать solve_bvp

from scipy.integrate import solve_bvp, solve_ivp
import matplotlib.pyplot as plt
import numpy as np

def eq(t,u): return [u[1], -0.106*u[1]-0.006*u[0]] #return u' and u''
def bc(u0,u1): return [u0[0], u1[1]-1 ]

res = solve_bvp(eq, bc, [0,1], [[0,1],[1,1]], tol=1e-3)
print res.message

# plot the piecewise polynomial interpolation, 
# using the last segment for extrapolation
time = np.linspace(0, 10, 501) 
plt.plot(time, res.sol(time)[0], '-', lw=2, label="BVP extrapolated")
# solve the extended solution as IVP solution on [1,10]
ivp = solve_ivp(eq, time[[0,-1]], res.y[:,0], t_eval=time)
plt.plot(time, ivp.y[0], '-', lw=2, label="IVP from BVP IC")

# plot decorations
plt.xticks(range(0,11))
plt.grid(True)
plt.legend()
plt.xlabel('time')
plt.ylabel('u(t)')
plt.show()

Обратите внимание, что продолжение выполняется путем экстраполяции из заданного интервала [0,1] до [0,10] и что значения в 1 имеют допуск 1e-3. Таким образом, можно получить лучший результат на большом интервале, используя solve_ivp с вычисленными значениями при t = 1 в качестве начальных значений. Разница в этом примере составляет около 0,01.