python: Initial condition in solving differential equation

Question

I want to solve this differential equation: y′′+2y′+2y=cos(2x) with initial conditions:

1. y(1)=2,y′(2)=0.5

2. y′(1)=1,y′(2)=0.8

3. y(1)=0,y(2)=1

and it's code is:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def dU_dx(U, x):
    return [U[1], -2*U[1] - 2*U[0] + np.cos(2*x)]
U0 = [1,0]
xs = np.linspace(0, 10, 200)
Us = odeint(dU_dx, U0, xs)
ys = Us[:,0]
plt.xlabel("x")
plt.ylabel("y")
plt.title("Damped harmonic oscillator")
plt.plot(xs,ys);

how can I fulfill it?

Answer 1

Your initial conditions are not, as they give values at two different points. These are all boundary conditions.

def bc1(u1,u2): return [u1[0]-2.0,u2[1]-0.5]
def bc2(u1,u2): return [u1[1]-1.0,u2[1]-0.8]
def bc3(u1,u2): return [u1[0]-0.0,u2[0]-1.0]

You need a BVP solver to solve these boundary value problems.

You can either make your own solver using the shooting method, in case 1 as

def shoot(b): return odeint(dU_dx,[2,b],[1,2])[-1,1]-0.5

b = fsolve(shoot,0) 

T = linspace(1,2,N)
U = odeint(dU_dx,[2,b],T)

or use the secant method instead of scipy.optimize.fsolve, as the problem is linear this should converge in 1, at most 2 steps.

Or you can use the scipy.integrate.solve_bvp solver (which is perhaps newer than the question?). Your task is similar to the documented examples. Note that the argument order in the ODE function is switched in all other solvers, even in odeint you can give the option tfirst=True.

def dudx(x,u): return [u[1], np.cos(2*x)-2*(u[1]+u[0])]

---

Solutions generated with solve_bvp, the nodes are the automatically generated subdivision of the integration interval, their density tells how "non-flat" the ODE is in that region.

xplot=np.linspace(1,2,161)
for k,bc in enumerate([bc1,bc2,bc3]):
    res = solve_bvp(dudx, bc, [1.0,2.0], [[0,0],[0,0]], tol=1e-5)
    print res.message
    l,=plt.plot(res.x,res.y[0],'x')
    c = l.get_color()
    plt.plot(xplot, res.sol(xplot)[0],c=c, label="%d."%(k+1))

---

Solutions generated using the shooting method using the initial values at x=0 as unknown parameters to then obtain the solution trajectories for the interval [0,3].

x = np.linspace(0,3,301)
for k,bc in enumerate([bc1,bc2,bc3]):
    def shoot(u0): u = odeint(dudx,u0,[0,1,2],tfirst=True); return bc(u[1],u[2])
    u0 = fsolve(shoot,[0,0])
    u = odeint(dudx,u0,x,tfirst=True);
    l, = plt.plot(x, u[:,0], label="%d."%(k+1))
    c = l.get_color()
    plt.plot(x[::100],u[::100,0],'x',c=c)