How to write initial conditions in scipy.integrate.ode (or another) function?

Question

I'm trying to solve differential equation using python scipy.integrate.odeint function and compare it to the mathcad solution.

So my equition is u'' + 0.106u'+ 0.006u = 0, the problem I'm stuck in is the initial conditions which are u(0)=0 and u'(1)=1. I don't understand how to set u'(1)=1.

Python code:

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 code:

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

Mathcad diagram looks like this:

https://pp.userapi.com/c850032/v850032634/108079/He1JsQonhpk.jpg
which is etalon.

And my python output is:

https://pp.userapi.com/c850032/v850032634/10809c/KB_HDekc8Fk.jpg
which does look similar, but clearly the u(t) is incorrect.

Answer 1

This is a boundary value problem, you need to use 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()

Note that the continuation is by extrapolation from the given interval [0,1] to [0,10] and that the values at 1 are with a tolerance of 1e-3. So one could get a better result over the large interval by using solve_ivp with the computed values at t=1 as initial values. The difference in this example is about 0.01.